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ABSTRACT 

Post-test  analysis  of  vehicle  tracking  data  taken  from  a  three  dimensional  underA  a- 
ter  range  requires  the  splicing  together  of  several  pieces  of  track  individually  produced 
by  distinct  transducer  arrays.  The  locations  and  orientations  of  the  arrays  must  be 
known  precisely  in  order  to  convert  locally  determined  tracks  into  a  coherent  record  in 
the  general  range  coordinate  system.  The  maintenance  of  calibration  of  the  system  is  a 
problem  and  this  study  presents  some  least  squares  methodology  that  uses  the  tracking 
data  itself  to  monitor  it.  More  specificially,  methodologies  arc  developed  to  improve 
upon  the  array  displacement  and  orientation  correction  algorithms  and  to  estimate  tim¬ 
ing  synchronization  offsets.  Applications  are  made  to  real  data. 
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I.  INTRODUCTION  AND  BACKGROUND 
A.  BACKGROUND 

Torpedoes  are  tested  for  proper  operation  at  the  Dabob  Bay  and  Nanoose  under¬ 
water  ranges  operated  by  the  Naval  Undersea  Warfare  Engineering  Station,  Keyport, 
Wa.  (NUWES).  There,  the  torpedo  is  tracked  by  a. system  of  liydrophone  arrays  on  the 
ocean  floor.  In  order  to  determine  if  tlie  torpedo  is  acquiring  and  homing  on  its  target 
properly,  accurate  tracking  is  essential. 

The  arrays  in  the  system  are  of  the  short  base  line  type,  each  containing  four  sonar 
transducers  placed  rigidly  at  the  corners  of  a  cube  in  a  manner  that  describes  a  cartesian 
coordinate  system  in  three  dimensions  (See  Figure  1).  An  array  receives  a  distinctive 
signal  from  a  synchronously  timed  pingcr  attached  to  the  torpedo.  The  differentials  of 
the  sound  wave  front  times  of  arrival  at  tlie  four  hydrophones  allow  the  coinjmtation 
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of  the  azimuth  and  elevation  angles  of  the  normal  to  the  wave  front  at  the  origin  of  the 
local  coordinate  system.  Then  assuming  direct  path  propagation,  one  can  ray  trace  us¬ 
ing  Snell's  law  (Ref  1).  starting  with  the  aforementioned  elevation  angle  and  utilizing  a 
velocity-versus-depth  profile  for  the  speed  of  sound  in  the  water.  Finally,  the  time  dif¬ 
ferential  between  the  source  pulse  at  the  torpedo  and  its  arrival  at  the  array  is  used  to 
stop  the  ray-tracing  algorithm  and  determine  the  location  of  the  torpedo  relative  to  the 
array.  The  local  track  is  the  sequential  set  of  these  estimated  positions. 

The  functioning  of  the  Range's  tracking  system  requires  knowledge  of  the  location 
and  orientation  of  each  array.  Moreover  once  these  values  are  established,  they  must 
be  monitored  against  slippages  of  various  forms,  i.e.,  calibration  must  be  maintained. 
Also,  the  synchronization  of  the  pinger  with  the  shore  based  computers  requires  great 
precision  .  because  timing  errors  can  be  confused  with  other  calibration  errors. 

B.  INTRODUCTION 

Each  array  in  the  range  system  operates  over  a  limited  radius.  As  the  target  passes 
through  the  range,  it  is  tracked  at  discrete  time  points  (called  point  counts)  by  a  number 
of  these  arrays  (See  Figure  2  on  page  3).  The  overall  path  is  constructed  by  trans¬ 
forming  each  piece  of  local  track  to  the  coordinates  of  the  range  based  upon  the  assumed 
location  and  orientation  of  the  various  local  coordinate  systems.  Discontinuities,  or 
mismatches  occur  because  the  track  produced  by  one  array  does  not  mesh  exactly  with 
that  produced  by  a  neighboring  array.  If  the  tracking  were  free  of  noise,  but  the  actual 
positions  or  orientations  of  the  arrays  were  dilTerent  from  the  ones  assumed,  errors  in 
the  two  versions  of  the  track  would  be  obvious  and  interpretable. 

A  second  possible  source  of  error  is  due  to  the  water  conditions.  The  ray-tracing 
procedure  depends  on  the  velocity  of  sound  in  each  layer  of  water.  A  measurement  of 
the  sound  velocity,  from  the  surface  to  the  bottom,  is  taken  for  each  day's  operation. 
This  information,  called  a  depth-velocity  profile,  is  assumed  to  be  constant  throughout 
the  range  for  that  day's  operations.  The  velocity  of  sound  in  water  is  highly  variable, 
and  influenced  by  depth,  temperature,  and  salinity  changes.  Any  errors,  either  a  bias 
of  the  measurement  or  inhomogeneity  from  one  part  of  the  range  to  another,  will  affect 
the  accuracy  of  the  range.  Also  any  noise  in  the  water  will  influence  the  accuracy  of  the 
range. 
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Figure  2.  Plan  view  of  range  coordinate  system  and  array  tracking  regions. 


The  purpose  of  this  thesis  is  to  study  the  use  of  crossover  data  (i.e.,  two  versions  of 
a  track  for  tlie  same  point  count  set  \vh:n  the  vclticlc  is  in  the  range  of  two  or  more 
diiTercnt  arrays)  for  monitoring  the  calibration  of  the  torpedo  test  range.  A  major  goal 
is  to  address  the  question  of  whether  mismatches  in  the  dilferent  array  tracks  can  be 
attributed  to  slippages  in  array  position  and  orientation,  timing  errors  or  whctlicr  some 
other  source  of  systematic  error  should  be  treated. 

C.  STRUCTURE 

The  structure  of  the  remainder  of  the  thesis  is  as  follows:  Chapter  2  investigates  the 
source  of  error  in  the  array  locations.  Chapter  3  deals  with  the  timing  olTsct  error  be¬ 
tween  the  torpedo  pinger  and  the  shore  computer.  Chapter  4  gives  a  mctiiodology  in 
order  to  combine  multiple  array  displacement  and  orientation  corrections  into  a  single 
correction.  Conclusions  and  rccomcndations  are  provided  in  chapter  5. 


3 


II.  ARRAY  DISPLACEMENT  AND  ORIENTATION  CORRECTIONS 
A.  MODEL 

Consider  a  set  of  point  counts  S  in  a  crossover  region.  It  is  convenient  to  refer  to 
the  two  arrays  as  the  “left"  array  and  the  "right"  array.  Let  the  crossover  data  in  range 
coordinates  be 


Y{t) 


>2(0 

>3(0 


for  track  determined  by  the  left  array 


and 


X(t) 


-*^2(0 

J^3(0 


for  track  determined  by  the  right  array  for  t  e  S 


T  =  Sumber  of  points  in  S 

provided  by  two  different  sensing  arrays.  Let  us  agree  that  the  Y(t)  data  comes  from  the 
array  whose  location  and  orientation  are  established,  and  our  goal  is  to  check  the  cali¬ 
bration  of  the  other  array.  In  particular,  does  there  exist  a  3  vector  A  and  a  3  x  3 
orthonormal  matrix  5  #  /  such  that  the  adjusted  track  values, 

.\V)  =  A  +  BXii)  (2.1) 

are  in  better  agreement  with  the  Y(t)  than  are  unmodified  X(t)  ? 

The  vector  A  is  related  to  a  displacement  of  the  sensing  array  and  the  matrix  B  is 
related  to  a  correction  of  its  orientation.  If  we  let  a  be  the  location  of  the  array  in  range 
coordinates,  ^{i)  be  the  local  track  determined  by  the  array  relative  to  its  own  position 
and  orientation,  and  /?  be  the  orthonormal  orientation  adjustment  which  rotates  the  lo¬ 
cal  coordinates  into  directional  alignment  with  the  range  coordinates,  then 

X{t)  =  a  -h  pm  (2-2) 

and 
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.V(/)  =  .-i  +  5a  +  5/?ct/)  (2.3) 

The  corrected  location  and  orientation  adjustments  are  A  +  Ba  and  Bp,  respectively. 

Estimates  of  A  and  B  are  obtainable  using  the  principle  of  least  squares.  We  will 
minimize  the  average  square  deviation  between  Y{t)  and  A  +  BX(t)  for  each  point  count. 
Using  the  squared  norm  notation,  define  an  objective  function 

2  =  aye  H  }-(0  -  II  ^ 

i.e., 


There  is  a  constraint  that  accompanies  the  minimization  of  Q.  Recall  that  the 
sensing  arrays  are  rigid  local  cartesian  coordinate  systems.  If  they  have  slipped  (i.e.  , 
moved  physically  such  as  by  the  action  of  a  ship's  anchor,  hooking  a  cable)  then  they 
undergo  a  displacement  and  a  reorientation.  The  matrix  B  should  be  orthonormal  as  is 
the  matrix  /?  in  the  original  calibration,  thus  our  problem  can  be  written  as. 

min  Q  =  aye  ||  T(/)  -  A  —  BX(i)  11^  (2.5) 

subject  to 

5^5  =  /  (2.6) 

where  I  is  a  3  x  3  identity  matrix. 

Solution  methodology  and  a  program  (Keymain)  for  this  procedure  are  presented  in 
[Ref.  2). 

B.  CROSSOVER  DATA 

The  crossover  data  described  above  forms  but  a  small  part  of  the  data  records  col¬ 
lected  during  regular  testing  operations.  Our  purpose  requires  software  to  extract  and 
organize  the  crossover  data  portions. 

The  Keyport  torpedo  tracking  data  file  (which  is  also  called  a  T  file)  is  a  list  of  re¬ 
cords.  Each  record  includes  point  count,  array  number  and  the  position  of  the  torpedo 
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in  the  three  dimensional  range  coordinate  system.  There  are  approximately  3000  records 
in  a  file.  Most  of  the  records  in  a  file  are  single  zone  records.  In  other  words  there  is 
only  one  array  that  receives  the  torpedo  signal  at  that  moment  (point  count). 

In  order  to  seperate  crossover  data,  that  is  the  group  of  records  taken  from  two  or 
more  different  arrays  at  the  same  time,  we  used  the  program  Keygate  which  was  devel¬ 
oped  under  the  direction  of  Prof  Robert  R.  Read,  Naval  Postgraduate  School,  1985. 
Fortran  77  code  for  this  program  is  presented  in  appendix  A.  The  21  August  1987 
Keyport  T  file  was  selected  for  illustration.  The  program  Keygate  was  used  to  produce 
a  set  of  crossover  data,  provided  by  the  arrays  3,  4  and  13.  Figure  3  on  page  7  shows 
such  a  set  and  provides  a  general  plan  view  in  a  macro  scale.  It  shows  array  locations 
and  torpedo  tracks  from  array  3,  array  4  and  array  13  in  the  crossover  zone.  Since  the 
coverage  of  the  plot  is  veiy  large,  it  is  too  difficult  to  discriminate  the  three  determi¬ 
nations  of  track.  This  particular  set  is  a  bit  rare  in  that  three  arrays  track  the  target 
vehicle  for  these  point  counts. 

Figure  4  on  page  8  gives  a  more  detailed  (micro)  view  for  these  tracks  in  the  cross¬ 
over  zone.  The  track  plotted  with  the  star  symbol  belongs  to  array  3.  The  other  two 
tracks,  which  are  plotted  with  plus  and  circle  symbols,  are  provided  by  array  4  and  array 
13  respectively.  Mismatches  between  the  tracks  can  be  easily  seen  in  this  figure. 

C.  TEMINOLOGY,  NOTATION  AND  RESULTING  CORRECTIONS 

Assuming  the  location  and  orientation  of  the  array  3  are  already  established,  we  can 
check  the  relative  calibration  of  array  4  and  array  13.  Recall  that  the  vector  A  is  related 
to  a  displacement  of  the  sensing  array  and  the  matrix  B  is  related  to  a  correction  of  its 
orientation.  The  estimated  displacement  vector  (D)  can  be  obtained  from  (2.3),  i.e., 

D  =  A  +  {B-r)0L  (2.7) 

To  describe  the  orientation  corrections  we  use  the  Euler  angle  representation  of  B  as 
developed  in  [Ref.  2,3).  Let  us  shorten  the  writing:  let  c,  =  cos  0,  and  s,  =  sin  4),  for 
i=  1,2,3,  (where  <i>,  is  the  rotation  angle  for  axis  i);  then  the  individual  planar  rotations 
can  be  writen  as  follows: 

10  0  C2  0  -52  Cj  -53  0 

p,  =  0  Cl  -5,  P2=  0  \  0  p3  =  S3  C3  0  (2.8) 

0  Si  Cj  S2  0  C2  0  0  1 
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For  example,  in  the  plane  of  (jc,  ,  jcj)  with  x,  held  fixed  the  effect  of  p,  is  the  usual  rota¬ 
tion  of  coordinates  (See  Figure  5  on  page  9J. 

These  three  rotation  angles  when  used  together  describe  the  three  degrees  of  freedom 
contained  in  the  matrix  B.  Specificially, 

5  =  P3P2Pi  (2.9) 

In  order  to  present  the  analysis  with  simplified  summaries,  the  three  components  of  the 
displacement  (2.7)  can  be  replaced  by  their  magnitude  (MD). 

MD  =  l  II D  II  (2.10) 
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Figure  4.  Plan  View  of  the  tracks  in  the  crossover  zone. 

Also  the  three  Euler  angles  can  be  replaced  by  one  maximum  angle  of  rotation  (See  Ref. 
3].  We  can  call  the  maximum  angle  of  rotation  the  total  rotation. 

Using  the  aforementioned  crossover  data  file  with  the  program  Keymain,  resulting 
summary  corrections  for  array  3  and  array  13  relative  to  array  2,  are  presented  in 
Table  1  on  page  9;  a  “worst  case"  has  been  chosen  deliberately.  In  this  table  QI  and 
QL  are  the  initial  and  the  last  values  of  the  least  squares  objective  function,  (2.4).  Also 
MD  and  TR  represent  the  magnitude  of  displacement  and  the  total  rotation  respectively. 

As  can  be  seen  from  Table  1  on  page  9,  corrections  for  the  array  locations  and  ori¬ 
entations  are  impossibly  large.  It  is  well  known  that  least  square  techniques  often  over 
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Figure  5.  Effect  of  rotation. 

optimize,  and  such  is  suspected  here.  Also  we  can  investigate  other  sources  of  error. 
We  will  deal  with  these  sources  in  the  following  chapters. 


Table  1.  INITIAL  DISPLACEMENT  AND  ORIENTATION  CORRECTIONS 


Array 

Qi  m 

QL  m 

%  Re¬ 
duction 

MD  (feet) 

TR  (de¬ 
gree) 

4 

18.32 

16.22 

11.04 

1577.94 

38.52 

in.  TIMING  OFFSET  ERROR 


A.  ANALYSIS 

Assume  there  is  a  timing  sjnchronization  error  between  the  torpedo  pinger  and  the 
shore  computer.  If  this  error  is  negative,  that  is  the  torpedo  pinger  sends  the  sound 
signal  prior  to  the  programmed  time,  then  the  calculated  distance  from  the  torpedo  will 
be  smaller  than  the  actual.  If  we  denote  the  timing  offset  error  by  A  and  the  speed  of 
sound  at  the  depth  of  the  torpedo  by  V'(t)  for  a  fixed  t  e  S,  the  distance  between  the  ac¬ 
tual  and  the  observed  position  of  the  target  will  be  AK(/).  If  the  timing  offset  error  is 
positive  then  the  calculated  distance  will  be  AI'(/)  larger  than  the  actual  ,  but  in  either 
case,  there  will  be  no  error  in  the  ray  trace  path  direction  in  the  vertical  plane  of  the 
array  and  the  torpedo. 

A  timing  offset  error  will  cause  the  track  determined  by  a  given  array  to  shift  (in  the 
horizontal)  along  the  direction  of  a  straight  line  from  the  array  to  the  target  for  each 
point  count.  To  illustrate,  without  cluttering  the  diagram,  four  point  counts  are  selected 
and  these  direction  lines  are  drawn  for  them  using  each  of  three  determinations  of  track 
by  the  three  arrays  in  Figure  6  on  page  11.  Note  that  if  the  target  determinations  are 
stretched  along  these  directions,  the  overall  coherence  of  the  three  versions  of  track  can 
be  improved.  There  will  also  be  some  stretching  in  the  vertical.  Although  it  is  nonlinear 
the  first  order  approximation  will  be  used. 

B.  METHODOLOGY 

1.  A  model  for  calculating  the  timing  offset  error 

Let  denote  the  observed  track  from  one  anay  in  the  range  coordinate  sys¬ 
tem. 

^""(0=  4(1)  (3.1) 


If  we  denote  the  array  location  in  the  range  coordinate  system  by  a 
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Figure  6.  Timing  offset  error  included  crossover  data 


(3.2) 


then  the  local  track  (with  orientation  the  same  as  that  of  the  range),  (A'^(r)),  for  this  ar¬ 
ray  can  be  calculated  by 


X\t) 


Jf{(0 

x[{t) 


=  X\i)-<x 


(3.3) 
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Let  be  the  azimuth  angle,  h{i)  be  the  horizontal  range  component  and  z{i) 
be  the  vertical  component  at  point  count  t.  Then  the  postion  of  the  torpedo  in  this  local 
coordinate  system  is 


/7(f)  cos  0(f) 

II 

4(0 

= 

h(r)  sin  0(f) 

4(0 

00 

If  there  is  a  timing  offset  error  of  size  A  (positive  or  negative),  then  the  correct 
local  track  will  be  shifted  (stretched  or  contracted)  along  the  ray  tracing  path.  A  first 
order  approximation  of  this  shift  requires  the  elevation  angle,  6{t),  of  the  ray  trace  in  the 
water  layer  of  the  observed  track.  Also  we  need  V'(t),  the  speed  of  sound  in  this  layer. 
Using  basic  trigonometn. ,  we  can  find  the  horizontal  and  vertical  components  of  the 
shift,  Ar(/)  cos  0(t)  and  AU(t)  sin  6(1)  respectively  (See  Figure  7  on  page  15).  Then  the 
corrected  horizontal  and  \ertical  components  are 

h{i,  A)  =  h{t)  +  AI'(r)  cos  d{i)  (3.5) 

z(t,  A)  =  r(t)  +  AUft)  sin  0(r)  (3.6) 


Using  Equation  (3.4),  we  can  write  the  corrected  position  in  the  local  coordi¬ 
nates. 


h{t)  cos  0(f) 

AU(f)  cos  0(f)  cos  0(f) 

h[i)  sin  0(f) 

AU(f)  cos  0(f)  sin  0(f) 

.  . 

AU(f)  sin  0(f) 

X^it,  A)  =  X^  + 


AV{t)  cos  6(t)  cos  0(f) 
AU(f)  cos  &({)  sin  0(f) 
A  U(f)  sin  0(t) 


(3.7) 


For  convenience,  let  a(t)  be  the  coefficient  of  the  error  term  in  the  local  coordinates, 


a(f)  = 


U(f)  cos  0(f)  cos  0(f) 
F(f)  cos  0(f)  sin  0(f) 
U(f)  sin  0(f) 


(3.8) 


so  we  can  write  the  track  in  local  coordinates  as  follows. 
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A'\i,  A)  =  X\t)  +  Aa(/) 

If  we  add  a  to  both  sides  of  the  above  equation,  then  the  track  in  range  corrdinates  can 
be  written  as  follows. 

X\t,  A)  =  X%)  +  Acj(/)  (3.9) 

where  is  the  corrected  track  in  the  range  coordinates.  If  we  use  the  notation 
A)  for  the  track  from  another  array, 

Y%,A)=-Y%)  +  Ab{i)  (3.10) 

where  b(t)  is  the  coefficient  of  the  error  term  for  the  second  array.  Corrected  tracks  in 
the  range  coordinate  system  from  both  arrays  are  given  by  (3.9)  and  (3.10). 

Now  we  can  calculate  the  timing  offset  error  by  minimizing  the  following  least 
squares  objective  function. 

minY  ll.V''(/.A)- P  (3.11) 

or 

min^  ll[A)  +  Aa(/)]-[y^0  +  A/>(/)]  11^  (3.12) 

After  taking  the  derivative  of  the  least  squares  objective  function  with  respect  to  A  and 
equating  to  zero,  we  can  find  an  explicit  solution  for  A. 

Y,^X%,  A)  -  y^(r,  A)] ^[<2(0  -  fc(0n 

A  = - -  (3.13) 

t 

In  order  to  see  the  result  is  a  minimum,  one  can  easily  take  the  second  derivative 
with  respect  to  A  and  check  that  it  is  larger  than  zero. 

For  calculating  a(t)  and  b(t)  we  have  to  find  V(t).  It  can  be  found  from  the 
depth-velocity  profile  .  Also  </>(/)  can  be  calculated  from  Eq.  (3.4) 
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4>{i)  =  tan' 


Xjlt) 

-V2(0 


(3.14) 


The  problem  here  is  to  find  the  elevation  angle  at  the  depth  of  torpedo  (0(/)). 
In  order  to  find  it,  one  has  to  use  a  ray  tracing  algorithm. 

2.  Ra>1racing 

a.  Background 

We  will  use  a  raytracing  algorithm  in  order  to  find  the  elevation  angle  at  the 
depth  of  torpedo.  There  are  two  general  models  used  in  raytracing.  Both  methods  in¬ 
volve  dividing  the  water  into  relatively  thin  layers,  and  then  modelling  the  ray  path 
through  each  succesive  layer. 

The  first  method  is  ISOSPEED.  In  this  model,  the  speed  of  sound  is  as¬ 
sumed  to  be  constant  in  each  layer.  The  velocity  chosen  is  the  mean  velocity  for  the 
layer.  The  ray  path  is  decribed  from  layer  to  layer  by  repeated  application  of  Snell's  law. 

The  second  method  is  ISOGRADIENT.  In  this  model  the  velocity  of 
sound  in  water  assumed  to  be  linear  and  described  by  the  follwing  equation. 


r,  =  1 0  -f  p,  X  z 


(3.15) 


where  I  \  is  the  speed  of  sound  at  the  water  surface,  Z  is  the  depth,  and  f',  is  the  gradi¬ 
ent.  The  gradient  represents  the  rate  of  change  in  velocity  as  the  depth  increases.  When 
using  a  model  with  a  single  layer,  and  P,  are  determined  by  least  squares  regression 
of  velocity  and  depth.  In  the  multilayer  model,  P,  is  the  velocity  extrapolated  to  the 
surface  (Z  <  0),  and  I',  is  the  gradient  within  that  layer. 

The  testing  of  both  approaches  proved  this  to  be  an  extremely  delicate  cal¬ 
culation.  An  algorithm  using  the  isospeed  method  was  less  analytically  complex  and 
computationally  less  intensive,  but  less  accurate  than  the  isogradient  case.  Therefore 
we  will  consider  only  isogradient  raytracing. 
b.  Isogradient  Raytracing 

To  begin  a  discussion  of  isogradient  raytracing,  assume  that  we  are  dealing 
with  only  one  layer,  and  that  the  velocity  of  sound  in  that  layer  is  given  by  Eq.  (3.I5). 
If  we  denote  the  position  of  the  array  and  the  position  of  the  torpedo  by  (/I, ,  Aj)  and 
(P^ ,  Pi)  in  the  vertical  plane  that  contains  the  array  and  the  torpedo,  then  the  ray  path 
is  a  circular  arc  in  that  plane  with  center  (C, ,  C,).  We  need  equations  to  find  the  center 
from  the  given  parameters.  Cj  is  given  by  the  following  equation. 
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Figure  7.  Sound  wave  path 


c  _ _ - 

Cl-  y 


(3.16) 


To  find  C„  let  R  be  the  length  of  the  radius  of  the  circle.  Calculating  the  distance  yields 
the  following  equations. 


{A^-C^f  +  {A2-C2f  =  R^ 


(3.17) 


(/>,-C,)'  +  (/>2-C2)'  =  /?' 


(3.18) 


Combining  these  two  equations  to  solve  for  C,  yields  the  following. 
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(3.19) 


C.= 


^1  +  ^1 
2 


2[P,~A,) 


+  -2Q 


Using  trigonometric  substitions,  C,  can  also  be  expressed  by  following  equation. 


(Cj  -  A^) 
tan  6 


(3.20) 


where  6  is  the  elevation  angle  at  the  array,  measured  from  the  horizontal. 

The  sound  velocity  does  not  remain  constant  along  this  arc.  Therefore,  the 
velocity  must  be  integrated  with  respect  to  distance  along  the  ray  path  to  determine  the 
time.  If  we  let  dg  be  the  elevation  angle  of  the  ray  at  the  torpedo,  and  be  the  elevation 
angle  of  the  ray  at  the  hydrophone,  then  the  time  of  transit  T  can  be  expressed  by  the 
following  equation. 

^  1  cos  0,(1  +  sin  0n) 


Other  quantities  needed  can  be  easily  obtained  by  trigonometric  substition 
into  these  equations  [See  Ref  4). 

c.  Multilayer  Raytracing 

In  actuality,  the  isogradient  approach  is  an  approximation  because  the 
gradient  is  not  linear.  Needed  accuracy  is  obtained  by  dividing  the  water  into  25  foot 
thick  layers.  Starting  at  the  ocean  surface  and  going  down  to  the  ocean  floor.  The  speed 
of  sound  is  measured  ever>'  25  feet.  The  difference  in  sound  velocity  from  one  meas¬ 
urement  to  another  is  used  to  determine  a  velocity  gradient  for  each  layer.  The 
isogradient  approach  is  then  used  within  each  layer.  Starting  with  a  guess  for  the  ele¬ 
vation  angle  at  the  array,  using  single  layer  approximation, 

01  =  tan  (3-22) 


we  can  calculate  the  ray  invariant  (RV).  Note  that,  the  ray  invariant  is  the  same  for  all 
layers. 


/?F  = 


cos  0, 

"vT 


(3.23) 
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where  is  the  speed  of  sound  at  the  depth  of  the  array. 

Let  I  be  the  layer  that  contains  tlie  target  and  J  be  the  layer  that  contains 
the  array.  If  we  arrange  layer  limits,  such  that  the  target  depth  will  be  the  upper  limit 
for  layer  I  and  the  depth  of  the  array  will  be  the  lower  limit  for  layer  J,  then  the  elevation 
angles  at  layer  I  and  at  layer  J  will  be  equal  to  the  elevation  angles  at  the  target  and  at 
the  array  respectively.  Note  that,  all  the  layers  between  I  and  J  have  the  thickness  of 
25  ft.  But  this  is  not  guaranteed  for  layer  I  and  layer  J. 

In  order  to  calculate  elevation  angle  at  the  depth  of  the  target  we  will  use 
a  recursive  algorithm.  First,  using  ray  invariant  we  can  calculate  the  elevation  angle  in 
each  layer  k  from  I  to  J. 

0j^  =  cos"'C/?F  X  Fj  (3.24) 

where  1  *  is  the  speed  of  sound  in  layer  K.  Then  using  single  layer  isogradient  approxi¬ 
mation,  we  can  calculate  the  horizontal  distance  (/?*)  in  each  layer,  beginning  from  layer 
J  through  layer  1.  Sunmiing  the  horizontal  distances  in  each  layer  from  1  to  J,  gives  the 
horizontal  distance  of  the  torpedo  from  the  array. 

y 

R  =  yRj  (3.25) 

k=l 

If  this  distance  is  in  the  tolerance  limits  of  the  observed  horizontal  range,  we  can  stop 
the  algorithm.  Otherwise  we  have  to  update  our  initial  guess  for  the  elevation  angle  at 
the  array  by  the  following  formula, 

=  (3.26) 

and  calculate  the  horizontal  range  until  reaching  the  tolorance  limits  of  observed  range. 
When  we  find  the  horizontal  range  from  the  algorithm  within  the  tolerance  limits  of  the 
observed  horizontal  range,  the  last  calculated  elevation  angle  in  layer  J  is  the  elevation 
angle  we  desired  [See  Ref  5]. 

3.  Program  TOE 

The  author's  Fortran  77  program  to  implement  the  procedure  for  finding  and 
eliminating  the  timing  offset  error,  is  presented  in  appendix  B.  It  uses  the  crossover  data 
from  the  output  of  program  Keygate.  Other  inputs  are  the  depth-velocity  profile  and 
the  range  configuration  file. 
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ui  Sr  \\ 

O  I 


o  .  ORIGINAL  TRACK  FROM  TRE  ARRAY  2 
•  -  CORRECTED  TRACK  FROM  THE  ARRAY  2 
0  -  ORIGINAL  TRACK  FROM  THE  ARRAY  1  ^ 
u  -  CORRECTED  TRACK  FROM  THE  ARRAY  1 1 


° . . ° . ° . . . . . . 0 . Q-...., 


24600 

DOWN  RANGE 


Figure  8.  Plan  view  for  the  original  and  corrected  tracks 


The  main  program  calls  five  subroutines.  Subroutine  Local  finds  the  local  track, 
than  subroutine  Transfrm  transforms  the  local  track  to  its  azimuth,  horizontal  and  ver¬ 
tical  components.  Subroutine  Lsline  calculates  the  speed  of  sound  at  the  water  surface 
and  the  gradient,  using  single  layer  approximation.  Subroutine  Rtrace  finds  the  ele¬ 
vation  angle  at  the  torpedo  using  the  multilayer  isogradient  raytracing  algorithm. 
Finally  subroutine  Cdelta  calculates  the  timing  offset  error  and  writes  the  pair  of  clean 
track  to  a  user  identified  output  file.  For  input  crossover  data  that  includes  30  point 
counts,  it  takes  approximately  20  seconds  to  run  this  program  on  an  IBM  XT. 
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4.  Result 

In  order  to  see  the  elTect  of  eliminating  timing  offset  error  we  produced  ten 
other  crossover  data  files  using  program  Kcygate.  These  data  files  come  from  three 
days.  We  thought  that,  if  there  is  a  timing  offset  error  in  the  data  then  it  must  be  the 
same  for  each  operation  (day)  and  each  array.  We  calculated  the  timing  offset  error  for 
each  of  the  ten  data  files.  Data  dates  and  summary  results  are  presented  in  Table  2  on 
page  20 

It  can  be  seen  from  Table  2  on  page  20  that  the  timing  offset  error  is  fixed  for 
each  day.  Figure  8  on  page  18  gives  a  plan  view  for  the  corrected  and  uncorrected 
crossover  data.  We  denote  uncorrected  data  with  the  file  N10S87.00I  and  corrected  data 
with  the  file  C10S87.001. 

We  calculated  the  array  displacement  and  orientation  corrections  from  ten  cor¬ 
rected  crossover  data  files.  Using  array  2  as  the  base  array,  we  got  corrections  for  the 
arrays  3  and  11.  A  summary  of  these  corrections  is  pre.sented  in  Table  3  on  page  23. 
As  before  QI  and  QL  are  the  initial  and  the  last  values  of  the  least  squares  objective 
function.  Also  .MD  and  TR  represent  the  total  displacement  and  total  rotation  respec¬ 
tively. 

It  can  be  seen  in  Table  3  on  page  23,  there  are  still  big  array  displacement  and 
orientation  corrections,  even  though  there  is  no  timing  offset  error  in  the  data.  We  got 
large  amounts  especially  in  the  displacement  corrections  for  the  vertical  axsis.  Since  we 
know  that  there  is  no  effect  in  the  water  that  moves  the  array  in  the  vertical  direction, 
we  decided  to  add  another  condition  to  our  minimization  problem,  in  order  to  hold  the 
array  at  the  bottom  of  the  ocean  when  we  get  the  displacement  and  orientation  cor¬ 
rections. 

C,  A  MODIFICATION  TO  THE  RANGE  CALIBRATION  METHODOLOGY 
1.  Model 

Recall  from  chapter  2,  our  goal  is  to  estimate  A  and  B  using  the  principle  of 
least  squares.  Our  objective  function  was  stated  in  Eq.  (2.4)  and  Eq.  (2.5).  Since  we 
decided  to  hold  the  array  at  the  bottom,  we  have  to  add  another  constraint  to  the  min¬ 
imization  problem,  i.e.,  A  =  0- 

Our  goal  is  to  estimate  A  and  B  while  minimizing  the  least  square  objective 
function  Q,  see  (2.4), 

(?  =  aye  ||  Y{i)- A- BX{t)  ||^  (3.27) 
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subject  to  the  constraints 
BB^=1 

7)3  =  0  (3.28) 

where  Z)'’  =  (D, ,  Dj ,  Dj),  see  (2.7). 

Table  2.  TIMING  OFFSET  ERROR  FOR  THE  CROSSOVER  DATA  FILES 


No 

File  Name 

Date 

Left 

Ar¬ 

ray 

Right 

Ar¬ 

ray 

Point 

Count 

Timing  Error 

n 

N05A87.001 

5  Aug. 
1987 

2 

3 

25 

-0.003 

2 

N05AS7.002 

5  Aug. 
1987 

2 

3 

25 

-0.003 

3 

N05A87.003 

5  Aug. 
1987 

2 

3 

25 

-0.003 

D 

N05A87.004 

5  Aug. 
1987  i 

2 

11 

20 

-0.003 

5 

N05A87.005 

5  Aug.  j 
1987 

2 

11 

! 

20 

-0.003 

1 

6 

N05A87.006 

5  Au2. 
1987 

2 

1 

11 

20 

-0.003 

B 

N05A87.00' 

5  Aug. 
1987 

2 

11 

20 

-0.003 

8 

N21A87.001 

21  Aug. 
1987 

1 

2 

3 

23 

-0.001 

9 

N2IA87.002 

21  Aug. 
1987 

2 

3 

26 

-0.001 

10 

N10S87.001 

10  Sep 
1987 

2 

11 

17 

-0.001 

Now  let  us  make  a  change  and  work  with  the  deviations. 

Y^t)=Y{t)-y  (3.29) 

X^t)  =  AXi)-x  (3.30) 

where  x  =  ave  X{i)  and  y  =  ave  Y{i) .  Then  we  can  write  the  objective  function  (3.27) 
as  follows. 
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Q  =  aye  ||  YJ^t)  +y- A-  +  JF)  P  (3.31) 

or 

e  =  aye  II  [  1'^/)  -  +  Cy  -  ^  V  (3.32) 

It  can  be  shown  (Pythagorean  theorem)  that,  (3.321  is  equal  to  the  following. 

Q  =  aye  ||  Y/t)  -  BXJ,t)  |P  +  \\y-A-Bx  W'^  (3.33) 

Also  from  (2.7)  we  A  =  D-{B—  IJa.  If  we  use  notation  Q,  for  orientation  and  Q, 
for  location, 

e^  =  aye  ||y^/)-BA’^/)  11^  (3.34) 

2/=  ILF-^-fi3Fp=  lly-a-D-fi(J-a)  11^  (3.35) 

then  we  can  write  the  objective  function  as  the  following. 

Q^Qo  +  Qi  (3-36) 

In  order  to  minimize  Q  we  should  minimize  Q,  and  Q,.  Now  let  the  covariance  matrices 
be 

C„  =  aye[A'J(f)J</0] 

C^^  =  aye[yJ'(r)r^r)] 

q..  =  aye[rJ(r)AVOJ  (3.37) 

since  B^B  =  I,  Q„  can  be  written  as  the  following  (See  Ref.  2  ]. 

Qo  =  lr{Cyy  -ICy^B  +  (3.38) 

where  tr  is  the  trace  operator.  Since  the  trace  is  a  linear  operator,  minimizing  Q,  is  the 
same  as  maximizing  the  following: 

max(C^,.,5^)  (3.39) 

B 
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A  solution  methodology  for  this  problem  is  presented  in  [Ref.  2].  After  some  algebra 
Q,  can  be  written  as  the  following. 

lly-a-B(J-a)  V -2D'^[y  -  a  -  B{x  -  a)'}  +  D'^ D  (3.40) 

Taking  the  derivative  with  respect  to  D  (See  Ref  3;  (eq.  1)]  and  equating  to  zero  provides 
an  estimator  for  D  (3.41)  which  minimizes  Q,. 

b  =  [S~0L-  B{x-a)']  (3.41) 


Since  Dj  must  equal  zero  (the  array  does  not  leave  the  bottom),  we  can  modify  (3.41)  to 
give  the  following. 

pji  -  «!  -  ^By{xj  -  a;)-j 

(3.42) 


Z)  = 


3 


T2  -  “2  -  j^ByiXj  -  (Xj) 


y=i 

0 


2.  Solution  Algorithm 

In  order  to  solve  the  problem  initialize  A  and  B  as  follows. 


B=I 

3 

fVl  -«!  -«;)-[ 

7=1 

D  =  _  _ 

>2  ~  *2  ~  ~ 

7=1 

0 

then  solve  for  B  by  maximizing  the  following  equation  in  the  usual  way  [Ref  2], 


(3.43) 


max  tr  [_Cy^B^'] 


(3.44) 


Next  calculate  D  with  (3.42)  and  continue  the  algorithm  by  calculating  another  B.  If  the 
difference  between  the  old  and  new  B  is  large,  then  turn  to  the  beginning  and  continue, 
i.e.,  calculate  another  D  and  from  this  D  calculate  another  B.  Stop  the  algorithm  when 
D  and  B  stabilize.  Since  the  calculation  of  B  while  maximizing  Eq.  (3.47)  requires  a 
similar  recursive  algorithm,  this  whole  procedure  is  analytically  ver>’  complex  and  com¬ 
putationally  veiy^  intensi'^e. 
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An  APL  code  to  implement  this  procedure  is  presented  in  appendix  C.  Corr  is 
the  main  function,  it  calculates  D,  QI  and  QL,  orientation  matrix  B  is  calculated  by  the 
functions  Twodim  and  Bmax. 

3.  Resulting  Corrections 

After  correcting  the  tracks  for  timing  offset  errors,  we  ran  the  program  Keymain 
again.  Then  the  array  displacement  and  orientation  corrections  from  ten  crossover  data 
files  are  presented  in  Table  4  on  page  25  As  before  QI  and  QL  are  the  initial  and  the 
last  values  of  the  objective  function.  Ajt,  Ay  and  Az  are  the  displacement  corrections 
(feet).  Also  Euler  angle  corrections  are  denoted  by  <f>^  and  (^3  (radians). 

One  can  easily  see  in  Table  4  on  page  25  that  the  location  and  orientation  cor¬ 
rections  are  quite  small  for  each  array.  This  certainly  improves  the  stuation  depicted  in 
Table  1  on  page  9.  Note  also  there  are  some  changes  from  one  observation  to  another 
for  the  same  array.  Since  all  changes  are  verj'  small,  we  decided  that,  these  corrections 
are  affected  by  some  kind  of  noise.  In  order  to  calculate  stable  array  corrections,  in 
other  words  to  calibrate  the  range,  we  need  to  remove  this  noise  from  the  array  dis¬ 
placement  and  orientation  corrections.  We  will  try  to  find  these  stable  corrections  in  the 
next  chapter. 


Table  3.  DISPLACEMENT  AND  ORIENTATION  CORRECTIONS  FOR  COR¬ 
RECTED  DATA  FILES 


File 

Array 

QI 

QL 

•’.b  Re¬ 
duction 

MD  (feet) 

TR  (de¬ 
gree) 

C05A  87.001 

3 

288 

22.6 

92.7 

32.9 

0.5 

C05A87.002 

169 

20.5 

87.8 

105.5 

1.6 

C05A87.003 

3 

51.8 

15.2 

70.6 

7.2 

0.2 

C05A87.004 

11 

36 

16.9 

53 

5503.7 

82.0 

C05A87.005 

11 

37.1 

12.8 

65.4 

1752.7 

24.1 

C05A87.006 

11 

53.8 

13.8 

74.3 

8390.7 

175.7 

C05A87.007 

11 

108 

17.5 

83.7 

5799.6 

83.7 

C21A87.001 

3 

166 

57.1 

65.6 

95.5 

1.06 

C21A87.002 

3 

295 

55.5 

81.1 

2314.2 

141.7 

C10S87.001  11  39.6  26.5  33  200.2 


But  before  going  on  let  us  gather  together  some  summaries  of  the  several  effects 
for  our  ten  cases.  These  are  provided  in  Table  5  on  page  25.  In  this  table  Raw  Ql  is 
the  average  difference  between  the  two  versions  of  uncorrected  track: 

Raw  QI  =  aye  It }'(/)  ~  ^(0  II  ^  (3-45) 

These  figures  have  the  dimension  of  square  feet  and  are  generally  quite  large.  The  root 
mean  square  values  (e.q.  .^1946.301  =44.1  feet  for  the  first  one)  provide  an  initial  raw 
figure  for  the  seperation  of  the  two  versions  of  track.  The  Raw  QIM  column  is  the  mean 
square  deviation  if  displacement  corrections  are  made  without  any  array  orientation 
changes.  These  values  provide  substantial  drops  from  the  Raw  QI  values.  If,  in  addition 
to  displacement  corrections,  we  also  make  orientation  corrections  then  further  im¬ 
provement  ensues,  but  not  much,  as  shown  by  the  Raw  QL  column. 

Next  let  us  see  what  can  be  accomplished  when  one  first  corrects  the  tracks  for 
timing  offset  errors.  The  resulting  values  of  (3.45)  appear  in  the  Cor  Ql  column.  (See 
also  QI  in  Table  3).  Clearly  these  values  represent  a  considerable  improvement  over  the 
Raw  QI  values.  When  displacement  and  orientation  corrections  (program  Keymain)  are 
applied  to  the  timing  corrected  data,  the  resulting  mean  square  deviation  values  appear 
in  the  Cor  QL  column.  These  values  compare  with  the  Raw  QL  values  and,  like  those 
values,  represent  displacement  and  orientation  corrections  that  often  are  impossibly 
large.  Finally  if  we  require  Dj  =  0,  i.e.,  the  array  can  not  leave  the  bottom,  the  final 
mean  square  deviations  are  marked  Cor  QF  (see  also  QL  in  Table  4).  The  important 
point  is  that  these  values  are  but  modestly  larger  than  the  Cor  QL  values,  and  at  the 
same  time  represent  but  small  array  corrections  as  evidenced  in  Table  4.  Thus  the 
combination  of  corrections  for  timing  errors  followed  by  array  corrections  that  do  not 
allow  the  array  to  leave  the  bottom  provide  us  with  a  very  tenable  set  of  values  for  our 
objective  function. 


Table  4.  FINAL  ARRAY  DISPLACEMENT  AND  ORIENTATION  COR¬ 
RECTIONS 


File  Name 

Ar¬ 

ray 

Qi 

QL 

Ax 

Ay 

Az 

4>i 

<f>2 

</>3 

C05A87,001 

3 

288 

24 

0.37 

0.02 

0 

-8.8E-3 

-6.7E-5 

2.7E-4 

C05A  87.002 

3 

169 

24 

0.31 

0.01 

0 

-9.6E-3 

-5.7E-5 

2.6E-4 

C05A87.003 

3 

51.8 

15.5 

0.07 

8.6E-3 

0 

5.1E-4 

-2.0E-4 

2.2E-4 

C05A87.004 

11 

36 

34.7 

1.5 

-0.1 

0 

-4.4E-4 

7.9E-7 

-2.2E-5 

C05A87.005 

11 

37.1 

14 

2.6 

-0.2 

0 

-2.0E-3 

-2.3E-6 

-2.1E-5 

C05A87.006 

11 

53.8 

17.2 

2.9 

-0.2 

0 

-2.5E-3 

-2.6E-6 

-1.6E-5 

C05A87.007 

11 

108 

20.5 

4.2 

-0.4 

0 

-3.9E-3 

-4.8E-7 

-9.7E-6 

C21AS7.001 

3 

166 

78 

0.1 

-0.02 

0 

0.01 

-3.7E-4 

3.1E-4 

C21A87.002 

3 

295 

107 

-0.4 

0.04 

0 

0.01 

-4.7E-4 

8.3E-5 

C10S87.001 

11 

39.6 

37.6 

1.4 

0.1 

0 

-5.8E-4 

-5.9E-6 

-2.2E-5 

Table  5.  COMPARISON  OF  THE  SEVERAL  EFFECTS 


File  Name 

Raw  QI 

Raw 

QIM 

Raw  QL 

Cor  QI 

Cor  QL 

Cor  QF 

05A87.00I 

1946.306 

26.753 

23.081 

288 

22.6 

24 

05A87.00: 

2062.929 

35.198 

20.289 

169 

20.5 

24 

05A87.003 

2197.208 

15.361 

15.313 

51.8 

15.2 

15.5 

05.A87.004 

483.882 

35.220 

16.510 

36 

16.9 

34.7 

05AS7.005 

548.032 

14.098 

12.607 

37.1 

12.8 

14 

05.A87.006 

630.1885 

17.429 

14.015 

53.8 

13.8 

17.2 

05AS7.()07 

779.075 

20.595 

17.415 

1 08 

17.5 

20.5 

21A87.0()1 

472.3^6 

74.759 

57.225 

166 

75.1 

78 

2IA87.002 

554.389 

103.397 

55.423 

295 

55.5 

107 

10S87.001 

183.488 

36.821 

26.428 

39.6 

26.5 

37,6 
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IV.  COMBINING  MULTIPLE  ARRAY  CORRECTIONS  INTO  ONE 

STABLE  CORRECTION 

A.  RANGE  CALIBRATION 

Final  displacement  and  orientation  corrections  for  array  3  and  array  II  are  pre¬ 
sented  in  Table  4  on  page  25.  It  can  be  seen  from  this  table  that  array  corrections  are 
quite  small  and  reasonable,  but  if  we  carefully  take  a  look  at  the  table,  we  can  see  some 
small  differences  between  the  corrections  for  the  same  array.  For  our  purposes  array 
displacement  and  orientation  corrections  must  be  fixed  for  each  array.  Since  there  is  no 
effect  in  the  water  to  change  the  array  location  and  orientations  in  a  couple  days,  there 
must  be  some  kind  of  noise  which  alTect  the  array  corrections.  If  this  is  true,  in  other 
words,  if  observed  corrections  include  additional  noise,  we  should  remove  this  noise  from 
the  corrections  in  order  to  find  stable  array  displacement  and  orientation  corrections. 
Removing  this  noise  leads  to  our  best  estimate  of  stable  corrections  for  each  array.  We 
can  calibrate  the  range  after  finding  stable  corrections  for  each  array  in  the  range.  Our 
approach  is  to  remove  the  noise  from  the  corrections  by  means  of  a  linear  mixed  model. 
We  will  deal  with  this  subject  in  the  following. 

B.  LINEAR  MODEL 
I.  Model 

Let  us  assume  that  the  observed  corrections  for  an  array  include  two  effects. 
One  of  them  is  the  stable  correction  for  that  array  and  the  other  one  is  a  day  effect  for 
each  day's  experimental  setup.  If  we  let  a,  be  the  stable  array  corrections  for  the  array 
i,  p,  be  the  elTect  for  day  j,  andy,^*  be  the  kth  set  of  observed  corrections  for  array  i  from 
day  j,  then  we  can  write 

y'ijk  =  “(  +  +  ^ijk  0 


where  are  the  random  measurement  errors, 

/  =  1,2,...,;? 
y=l,2,...,9 

(4.2) 
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where  is  the  number  of  observations  for  array  i  and  for  day  j.  Since  we  assume  that, 
the  array  effect  (fi)  is  random,  then  our  model  is  a  mixed  effect  model. 

Since  our  model  is  a  mixed  effects  model  and  our  data  is  unbalanced,  i.e.,  we 
do  not  have  equal  number  of  observations  for  each  day  and  for  each  array,  an  acceptable 
method  for  estimating  the  variance  components  of  this  model  is  the  fitting  constants 
method.  It  is  also  known  as  Henderson's  method  III  [See  Ref  6].  But  first  we  must 
recognize  that  our  corrections  are  multivariate. 

2.  Multivariate  Analysis 

Since  each  observation  is  a  six  element  vector,  consisting  of  six  corrections  for 
the  array,  our  analysis  will  be  a  multivariate  analysis.  There  will  be  two  covariance 
matrices  to  be  estimated,  namely,  the  between-group  and  the  within-group  covariance 
matrices.  These  two  covariance  matrices  are  most  often  estimated  by  forming  a 
multivariate  analysis  of  variance  and  equating  mean  square  matrices  to  their  expecta¬ 
tions. 

Now  let  us  turn  to  our  problem  and  write  our  model  again. 

}ijk  —  +  l^j  +  ^ijk 

i  =  l,2....,p 
7=  12,. ..,q 

(4.3) 

where  y,,*  is  a  6  x  1  vector  of  observations,  a,  is  a  6  x  I  vector  of  stable  corrections  for 
array  i,  is  a  6  x  1  vector  of  noise  for  day  j  and  is  a  6  x  1  vector  of  standard  errors. 
We  assume  that  the  are  independently  and  identically  distributed  (0,  L^)  random  vec¬ 
tors  which  are  independent  of  e,j^.  The  {«„»}  are  assumed  to  be  independently  and  iden¬ 
tically  distributed  (0,  L,)  random  vectors.  Using  a  modified  fitting-of-constants  method 
[See  Ref  7j  we  can  obtain  estimators  for  the  between-group  covariance  matrix  (MS^)  and 
the  with-in-group  covariance  matrix  (iV/S„),  which  are  approximately  unbiased, 

ip  =  c-\MS,-MS^)  (4.4) 

i,=  MS^  (4.5) 

where 
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(4.6) 


If  we  let 


p  q  p  q  fij 


y-  =  (Z  ZV)~'Z  Z  Z/^* 

^/=I  /=I  /  fel  /=!  A:=l 


(4.7) 


and 


(4.8) 


then  the  estimators  for  the  between-group  covariance  matrix  (  A/Sj  )  and  the  with-in- 
group  covariance  matrix  A/S.  are  the  following 


p 

MS,  =  {q  -1)-'^CVV',  -A"  .)tv/  ~y.f 
J=\ 


(4.9) 


-  C(?-1)  +  (P-1)]|"'V  H-IO) 

^'=>  /=1  -  'fcl'  fc'  fcf 

Both  matrices  from  these  estimators  are  nonnegative  definite,  but  this  is  not  guaranteed 
for  the  dilTerence.  There  is  a  high  probability  that  the  result  of  Eq.  (4.4)  is  not  nonneg¬ 
ative  definite  and  it  is  not  a  proper  covariance  matrix.  An  estimation  methodology 
which  produces  nonnegative  definite  covariance  matrices  is  presented  in  Ref.  [8]. 

If  MSi  —  A/S.  is  not  nonnegative  definite,  then  in  order  to  get  a  nonnegative 
definite  estimator  for  the  covariance  matrix  let  /,  ^  >  ...  >  A.,  be  the  ordered  roots 

of  the  following  equation. 

|A/Sj-/!iV/S^|  =0  (4.11) 


The  roots  A,  are  called  the  characteristic  roots  of  A/S,  in  the  metric  of  A/S*.  In  order  to 
find  /„  let  us  define  a  matrix  L  satisfying  following  equation. 


L^MS^L  =  I  (4.12) 

We  can  find  a  matrix  L  satisfying  (4.12)  by  Cholesky  decomposition.  Another  choice 
for  L  satisfying  (4.12)  is 

L  =  =  Q^diag  {vr°'  ...  qI  (4.13) 

where  v,  are  the  eigenvalues  of  MS,  and  Q,  is  the  matrix  of  the  corresponding 
orthonormal  eigenvectors  of  .1/S^.  With  any  choice  of  L  satisfying  (4.12),  the  roots 
are  the  eigenvalues  of  a  symmetric  matrix  UMS^L  \ow  define  a  matrix  P  by 

P  =  L~^Q  (4.14) 

where  Q  is  the  matrix  of  the  corresponding  orthonormal  eigenvectors  of  UMSi,L  . 

If  >  1  for  all  i=  1,2,. ..,6  then  (A/S*  —  MS,)  is  nonnegative  definite  and  the  es¬ 
timator  (4.4)  is  in  the  parameter  space.  If  /,  <  1  for  some  i,  then  (.V/S*  —  MS,)  is  not 

nonnegative  definite.  Let  k  be  the  number  of  such  that  /.,>  1.  Furthermore  let 

A*  =  A;  =  diag  and  ?  =  (?*,  P,),  where  P*  is  6xk. 

Then  we  can  write, 

A/S,  -  A/S,  =  P,(A,  -  I,,)PI+  P^\,-I„)PI  (4.15) 

A/S,  -  A/S,  =  (A/S,  -  A/S,)+  +  (A/S,  -  A/S,).  (4.16) 

Hence  by  using  the  metric  of  A/S„,  we  have  written  (A/S*  —  MS,)  as  the  sum  of  a  non¬ 
negative  definite  matrix  (A/S*  —  MS,),  and  a  negative  definite  matrix  (A/S*  —  MS,)..  The 
(.V/S*  —  A/S„).  can  be  considered  to  be  the  closest  matrix  to  (A/S*  — A/SJ  among  all 
nonnegative  definite  matrices  that  have  the  same  characteristic  vectors  in  the  metric  of 
MS,  as  those  of  (A/S*  -  .V/SJ.  In  this  sense  we  can  interperet  the  term  (A/S*  -  A/S„). 
as  the  nonnegative  definite  portion  of  (A/S*  —  MS,)  obtained  by  projecting 
(A/S*  -  A/S„)  in  the  metric  of  MS,  on  to  the  set  of  all  nonnegative  definite  matrices. 
Thus  it  is  natural  to  consider  (^V/S*  —  MS,),  as  an  estimator  of  cZ,  =  £C-V/S*  —  A/S.]  and 
the  remaining  part  of  A/S* 

.3/S,-(.V/S,-A/S,)+  (4.17) 

as  an  estimator  of  S,.  An  intuitive  estimator  of  Z,  is  a  weighted  average  of  the  quantity 
in  (4.17)  and  MS, 
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Based  on  the  foregoing  discussion,  we  define  estimators  of  Lg  and  1,  in  model 

(4.1)  by 

=  0  (/•/:  =  0 
=  if  \<k<6 

and 


2^.  =  C(9  -  1)  +  2]"’  {{q  -  DIMS,  -dpi  +  zA/SJ 
where 

-  Up-l)  +  {q-l)l  (4.18) 

The  author's  APL  code  to  implement  this  procedure  are  presented  in  appendix 
D.  This  program  calls  the  APL  eigenvalue  librar>’  function  Eigenr  in  order  to  calculate 
eigenvalues  and  eigenvectors  of  a  real  matrix. 

3.  Multivariate  Estimation. 

The  estimation  of  the  fixed  effects  requires  values  for  Lp  and  L,.  We  will  use  the 
esiimators  developed  in  the  previous  section.  The  least  squares  estimator  for  a  requires 
the  variance  of  the  observations.  To  develop  this  we  change  notation.  Let  N  be  the 
total  number  of  experiments  (i.e.,  vector  observations),  then  the  model  can  be  written 
as 

Y  =  A  a  +  Z/?  +  e 
where 

r=[:vy]  /=i,...,a;=i . 6 

Af  =  Ixikl  A  k=l,...,q 

Z=[zJ  /=1,...,A  /=!,...,/» 

^  =  le,l 


2  = 


r  P  t  n 

‘-i=l  ;=i  -■ 


/  =  I,..., A  j=  1,...,6 


2  =  [ot^-yD  k=\....,p 
^  /=1 . q  j=\,...fi 

P  Q 

i=l  7=1 

jf,*  is  one  if  the  ith  observation  provides  a  correction  vector  for  the  ith  array,  zero  oth¬ 
erwise  and  2,1  is  one  if  the  ith  observation  was  collected  on  the  1th  day,  zero  otherwise. 
Note  that  for  our  data  \  =  10,  p=  2,  q  =  3.  The  /?,  =  (^;i, ••■,/?*)  are  independent  random 

vectors  for  1=  1 . q  with  mean  zero  and  covariance  matrix  Zp.  The  e,  =  {e,, . ej  are 

independent  error  random  vectors  for  i=  1,...,N  with  mean  zero  and  covariance  matrix 
1,.  All  are  independent  of  all  <?,.  Clearly 

El  T]  =  Aa  +  Z£[^]  +  =  A’a  (4. 19) 

To  obtain  covariances  it  is  convenient  to  change  our  mental  attitude  and  think 
of  the  Y,  ft,  /?,  e  arrays  as  vectors  for  which  the  six  values  of  the  second  subscript  are 
repeated  after  each  advancement  in  the  first  subscript.  With  this  attitude,  the  model  can 
be  described  using  matrix  direct  product  notation.  In  what  follows  liberal  use  is  made 
of  the  direct  product  calculus  developed  in  [Ref  9,  sec.  8.8). 

}'  =  (A'®4)(40ft)  +  (Z®4)(/^®/?)  +  c  (4.20) 

where  /,  and  4  are  identity  matrices  of  order  p,  q  and  6  respectively.  Now  the 
covariance  matrix  V  of  Y  can  be  developed  as 

r  =  £[  y  -  ( A'®  4)a  ]  [  }•  -  (T®  4)ft  ]  ^ 

£=  £C(Z®4)/?  +  c][(Z®4)/?  +  ey 

V=  £C(Z®4)/?/?^(Z^®4)3  -(-  Elee^l 

(Z®4)(4®I^)(Z^®4)  (7^®!,) 

£=ZZ^®S^-f /^®S,  (4.21) 

where  /v  is  the  identity  matrix  of  order  N. 
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The  least  squares  estimator  for  a  is  the  Aitken  estimator 


a  =  [(A'^04)  Y  (4.22) 

It  is  useful  to  use  the  structure  of  the  experimental  design  in  order  to  simplify  the  com¬ 
putation.  The  components  of  A'  =  are  one  or  zero  according  to  whether  the  ith 
experiment  involves  the  kth  sonar  array,  k=  l,..,p.  The  value  of  Z  =  Cz,/]  is  one  or  zero 
according  to  whether  the  ith  experiment  occurs  on  the  1th  day,  1=  l,...,q.  Note  that 
Z^Z  is  a  q  by  q  diagonal  matrix  (n„...,M,)  where  n,  is  the  number  of  experiments  included 
for  the  1th  day.  The  matrix  ZZ^  is  a  block  diagonal  matrix  (4, .•••,•/„,).  where  is  a 
square  matrix  of  order  n,  and  of  its  components  are  one.  It  follows  that  V  is  block  di¬ 
agonal  with  blocks  possessing  the  general  form 

1 ;  =  40-e  + 

Because  of  this  structure.  V-'  will  have  a  like  structure  and  we  need  the  inverse  of  indi¬ 
vidual  blocks.  This  can  be  accomplished  by  assuming  the  form 

I-' =  I„®a  +  (4.23) 

setting 

[40S^  -f-  y„0I^][40't  +  =  404 

and  solving.  Thus 

(40S,)(40a)  (7„0I^)(40a)  +  +  {J„®I,p){J„®d„) 

=  I„®o,a  -f-  J„®lL^a  +  =  404 

The  solution  is 


(4.24) 

(4.25) 

Now,  the  Aitken  estimator  for  a  can  be  managed  through  partitioned  matrices. 
The  rows  of  X  are  partitioned  according  to  Call  a  typical  one  X„.  Then,  identi¬ 

fying  AA  from  (4.22), 
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AA  =  =  C  ....(a;"'0«)  +  3 


(4,26) 


and  the  structure  of  Ajy,  is  that  of  a  p  by  n  matrix  whose  kth  row  is  constant  and  equal 
to  the  number  of  times  (experiments)  the  kth  sonar  array  was  included  on  that  particular 
day. 

We  also  need  from  (4.22),  BB  =  whose  value  can  be  obtained 

from  (4.26)  by  summing  the  quantities  of  the  form  [(A70«)  +  iXlJ^)®d„']{X„®l^)  .  The 
result  is 


DB  =  q 


(4,27) 


which  is  a  square  matrix  of  order  6p.  It  can  be  inverted  numerically  provided  p  is  not 
very  large.  Then  the  Aitken  estimator  for  a  can  be  written  as  follows. 


a=-{BB)~\AA)Y 


(4.28) 


An  APL  code  to  implement  this  procedure  for  our  data  is  presented  in  appendix 
E.  This  program  calls  the  APL  library  function  Penrose  in  order  to  find  a  generalized 
inverse  of  a  matrix.  Using  this  program  and  previous  data,  we  calculated  the  stable 
displacement  and  orientation  corrections  for  the  array  3  and  for  the  array  1 1  relative  to 
the  array  2.  Resulting  displacement  and  orientation  corrections  are  presented  in 
Table  6. 


Table  6.  STABLE  CORRECTIONS  FOR  THE  ARRAYS  RELATIVE  TO  THE 
BASE  ARRAY. 


.Array 

A  X 

A  \ 

AZ 

01 

0: 

03 

3 

0.54S81 

-0.0500 

0 

-4.3E-4 

-1.9E-5 

3.30E-5 

1! 

0.11021 

-0.0121 

0 

0.01400 

-9.6E-5 

5.03E-5 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 

A.  CONCLUSIONS 

As  we  pointed  out  at  the  beginning,  the  purpose  of  this  thesis  is  to  decide  whether 
or  not  crossover  data  can  be  used  to  monitor  the  calibration  of  the  torpedo  test  range. 
We  used  data  from  the  three  arrays,  using  one  array  as  the  base  array  and  obtained 
corrections  for  the  other  two  arrays.  Since  the  range  consists  of  24  arrays,  the  calibrated 
part  of  the  range  is  but  a  small  portion. 

Although  the  calibrated  part  of  the  range  is  a  tiny  portion,  corrections  for  the  two 
arrays  gave  us  an  idea  for  the  whole  range.  Since  observed  corrections  for  the  third  and 
eleventh  array  from  the  uncorrected  crossover  data,  were  enormously  big.  we  began  to 
investigate  other  sources  of  error.  First  wc  calculated  the  timing  offset  error  by  means 
of  Least  Squares.  Although  eliminating  the  timing  offset  error  did  not  reduce  the 
amount  of  corrections  as  much  as  expected,  we  decided  that  it  is  creditable,  because  no 
matter  what  array  it  is.  timing  offset  error  was  the  same  for  each  operation  on  the  same 
day.  Also  based  upon  plots,  it  is  very  encouraging  to  see  that  the  elimination  of  timing 
errors  reduces  the  mismatches  between  tracks  from  difl'erent  arrays.  This  gives  an  op¬ 
portunity  to  correct  not  only  the  crossover  data  but  the  entire  torpedo  track. 

Next  we  tried  to  hold  the  array  at  the  bottom  of  the  ocean.  This  reduced  the 
amount  of  the  corrections,  also  it  seems  to  stabilize  the  result.  After  forcing  =  0, 
corrections  were  quite  small  but  not  fixed  for  each  array.  In  order  to  get  the  fixed  cor¬ 
rections,  we  decided  to  remove  the  noise  which  affected  the  corrections  by  means  of  a 
linear  mixed  model.  We  used  a  fitting  of  constants  method  in  order  to  estimate  the  fixed 
corrections.  Since  the  resulting  corrections  were  quite  small,  we  decided  that  the  known 
locations  of  the  two  arrays  were  correct  and  the  source  of  the  mismaches  is  the  timing 
offset  error  between  the  torpedo  and  the  shore  computer. 

B.  RECOMMENDATIONS 

In  this  study  we  got  corrections  for  two  arrays.  In  order  to  calibrate  the  whole 
range,  this  methodology  must  be  used  with  large  data  files.  This  requires  large  computer 
storage  and  long  CPU  time.  If  the  results  for  the  whole  range  are  as  good  as  the  results 
for  this  small  portion,  this  methodology  can  be  gathered  in  one  computer  program.  In 
order  to  do  this  some  decision  rules  are  required  such  as  what  amount  of  error  will  be 
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assigned  to  timing  ofiset  error  and  what  amount  of  error  will  be  assigned  to  calibration 
error. 

If  :he  ocean  bottom  is  flat,  holding  the  array  at  the  bottom  is  reasonable.  If  not, 
extra  work  must  be  done.  Because  of  this,  before  using  our  methodology  mapping  the 
ocean  bottom  is  required.  Also  this  methodology  uses  the  same  base  array  for  each 
operation.  A  modification  which  makes  different  base  arrays  acceptable  will  be  useful 
and  maximum  likelihood  estimation  can  be  used  to  estimate  fixed  array  corrections.  A 
multivariate  version  of  maximum  likelihood  estimation  is  presented  in  [Ref  10]  and  [Ref 
111- 
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APPENDIX  A.  PROGRAM  KEYGATE 


PROGRAM  KEYGATE 


***  Program  to  read  in  raw  data  from  Keyport  hydrophone  *** 
***  arrays,  segregate  it  by  mode,  and  throw  out  unusable  *** 
***  records.  The  output  of  this  program  is  to  be  read  in  *** 
***  by  the  program  KEYMAIN.  *** 


INTEGER*4  CRT,  KBD 
CHARACTER*13  DSNAME,  SITNAM 
CHARACTER*9  TEMPI,  TEMP2,  TEMP3,  TEMP4 


PARAMETER  (KBD=5,  CRT=6) 


WRITE(CRT,*)  '  Please  enter  the  name  of  your  input  file:  ' 
READCKBD, '(A)' )  DSNAME 

WRITE(CRT,*)  '  Please  enter  the  name  of  the  range' , 

+  '  configuration  file:  ' 

READ(KBD,'(A)')  SITNAM 


TEMPI  =  'TEMPI.  TMP’ 
TEMP2  =  'TEMP2.TMP' 
TEMP3  =  'TEMP3.TMP' 
TEMP4  =  'TEMP4.TMP' 


CALL  STRMODC CRT, KBD, TEMPI, DSNAME) 


CALL  PAIR(CRT,KBD,TEMP1,TEMP2) 


CALL  RANGE ( CRT , KBD , S ITNAM , TEMP  2 , TEMP3 ) 


CALL  PAIR2(CRT,KBD,TEMP3,TEMP4) 


CALL  REC(CRT,KBD,TEMP4) 


WRITE(CRT,*) '  Operation  complete.  KEYGATE  terminating...  ' 
C 

STOP 

END 

C 


36 


n n n o o o  oono  no  n  n  on  o  noooo  on 


SUBROUTINE  STRMOD( CRT, KBD, TEMP l.DSNAME)  _  _ 

**  PROGRAM  TO  STRIP  ALL  MODES  EXCEPT  2'S  AND  7'S  FROM  KEYPORT  ** 
**  DATA.  (2  INDICATES  TARGET  SHIP,  7  INDICATES  TORPEDO)  ** 

****VWf  ********************************************************** 


CHARACTER  D0*2,  DSNAME*13,  TEMP1*9 
INTEGER  PC,  ARRAY,  CRT,  NHEAD 

OPEN( 1 ,FILE=DSNAME, STATUS=' OLD* ) 


WRITE(CRT,*) '  How  many  records  of  header  do  you  want  to', 
+  '  strip  off  the  file?' 

RE AD (KBD,*)  NHEAD 
DO  11  I  =  1, NHEAD 
READd,*) 

11  CONTINUE 

WRITE(CRT,*) '  Input  mode  to  be  kept?' 

READCKBD,*)  NUM 

0PEN( 2 , F I LE=TEMP 1 , STATUS= ' NEW ' ) 


10  READ( 1,100, END=50 , ERR=40 ) PC , DO ,X , Y , Z , ARRAY , MODE 

IF(D0  .NE.  '  ')G0T0  20 

IFCMODE  .NE.  NUM)  GOTO  20 
WRITE( 2 , 1 10 ) PC ,X , Y , Z , ARRAY , MODE 
20  CONTINUE 

GOTO  10 

40  WRITE(CRT,*)'  THERE  WAS  AN  ERROR  IN  THE  FILE' 

50  CONTINUE 

100  F0RMAT(I5,A2,1X,F7.  1,2X,F7.  1,2X,F7.  1,30X,I2,2X,I1) 

110  F0RMAT(1X,I5,2X,3F10.  1,2X,I2,2X,I2) 

CLOSE  (UNIT=1) 

CLOSE (UNIT=2) 

RETURN 

END 


SUBROUTINE  PAIR(CRT,KBD,TEMP1,TEMP2) 


***  PROGRAM  TO  PAIR  POINT  COUNTS  AFTER  THE  DATA  HAS  BEEN  *** 
***  GATED  BY  "  STRMOD  SECOND  PASS.  *** 


DIMENSION  X(200),  Y(200),  Z(200) 

INTEGER*4  PC(200),  ARRAY(200),  MODE(200),  CRT,  HOLD 
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CHARACTER*9  TEMPI,  TEMP2 


OPEN( 1 , F I LE=TEMP 1 , STATUS= ' OLD ' ) 
OPEN( 2 ,FILE=TEMP2 , STATUS=’ NEW' ) 

HOLD  =  0 
IREC  =  0 
NREC  =  0 
1  =  1 


.  .  .  READ  RECORDS  BY  TWO’S  TO  COMPARE  POINT  COUNTS: 

READ(l,*,ENti=40,ERR=30)  PC(I),  X(I),  Y(I),  Z(I), 

+  ARRAY(I),  MODE(I) 

10  READ(1,*,END=40,ERR=30)  PC(I+1),  X(I+1).  Y(I+1),  Z(I+1), 

+  ARRAY(I+1),  M0DE(I+1) 

NREC  =  NREC  +  1 

IF(PC(I+1)  .EQ.  HOLD)  THEN 

WRITE(2,100)  PC(I+1),  X(I+1),  Y(I+1),  Z(I+1), 

+  ARRAY(I+1),  M0DE(I+1) 

HOLD  =  PC(I+1) 

GO  TO  20 
END  IF 

IF(PC(I)  .EQ.  PC(I+1))  THEN 

WRITE(2,100)  PC(I),  X(I),  Y(I),  Z(I),  ARRAY(I),  MODE(I) 
WRITE(2,100)  PC(I+1),  X(I+1),  Y(I+1),  Z(I+1), 

+  ARRAY(I+1),  M0DE(I+1) 

HOLD  =  PC(I+1) 

IREC  =  IREC  +  1 
GO  TO  20 
END  IF 

IF(PC(I)  .NE.  PC(I+1))  THEN 
PC(1)  =  PC(I+1) 

X(l)  =  X(I+1) 

Y(l)  =  Y(I+1) 

Z(l)  =  Z(I+1) 

ARRAY(l)  =  ARRAY(I+1) 

MODE(l)  =  M0DE(I+1) 

1  =  1 
END  IF 

GOTO  10 

20  1  =  1  +  1 
GO  TO  10 


30  WRITECCRT,*)'  THERE  IS  AN  ERROR  IN  THE  DATA  FILE  IN  REC. '.NREC 
WRITECCRT,*)'  ...  OPERATION  TERMINATING  DUE  TO  ERROR.  ' 

STOP 
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40  CONTINUE 

CL0SE(UNIT=1) 

CL0SE(UNIT=2) 

100  F0RMAT(1X,I5.2X,F7.  1,2X,2F9.  1,2X,I2,2X,I2) 

RETURN 

END 


SUBROUTINE  RANGE( CRT , KBD , SITNAM , TEMP2 , TEMP3 ) 


***  This  program  completes  the  third  gating  of  Keyport  *** 
***  range  data.  It  reads  array  location  data  from  a  *** 
*Vf*  site  specific  configuration  file  and  tests  to  see  if  the*** 
***  data  is  in  the  valid  overlap  area.  *** 
************■>'?********************★*********•********************* 


INTEGER*4  ARRAY,  CARRAY,  CRT,  PC 
REAL*4  CONI'IG(200,4),  LX,  LY,  LZ,  MAXVAL 
CHARACTER  SITNAM*13,  TEMP2*9,  TEMP3*9 


Open  input  and  output  files: 
0PEN( 1 ,FILE=TEMP2 , STATUS* 'OLD' ) 
0PEN(2,FILE=SITNAM,STATUS=’0LD’ ) 
OPEN(3,FILE=TEMP3,STATUS='NEW' ) 


.  .  .  Read  site  configuration  into  CONFIG  array: 

NREC  =  0 
1  =  1 

10  READ(2,*,END=40,ERR=30)  C0NFIG(I,1),  C0NFIG(I,2),  C0NFIG(I,3), 

+  C0NFIG(I,4) 

NREC  =  NREC  +  1 
1  =  1  +  1 
GOTO  10 

30  WRITE(CRT,*) '  There  was  an  error  reading  the  config' , 

+  '  file  in  record  ' ,NREC 

40  CONTINUE 

CL0SE(UNIT=2) 


NDREC  =  0 

...  Read  X,  Y,  Z,  and  ARRAY  from  input  data  file: 
45  READ(1,*,ERR=70,END=80)  PC,  X,  Y,  Z,  ARRAY,  MODE 

DO  50  I  =  l.NREC 
LX  =  CONFIGd,!) 
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LY  =  C0NFIG(I,2) 

LZ  =  C0NFIG(I,3) 

CARRAY  =  INT(C0NFIG(I,4)) 

MAXVAL  =  4700. 

. . .  Match  array  number  in  data  file  with  that  in  config.  file. 
...  If  they  are  equal  compute  slant  range  distance  (SRDIST): 
IF( ARRAY  .EQ.  CARRAY)  THEN 

SRDIST  =  SQRT(  (X  -  LX)**2  +  (Y  -  LY)**2  +  (Z  -  LZ)**2  ) 
IF(MAXVAL  .GE.  SRDIST)  THEN 
WRITE(3,100)  PC,  X,  y,  Z,  ARRAY,  MODE,  SRDIST 
NDREC  =  NDREC  +  1 
END  IF 
END  IF 

50  CONTINUE 

GOTO  45 

70  WRITE(CRT,*) '  There  is  an  error  in  the  data  file.  ' 

80  CONTINUE 

CLOSE ( UN IT=1) 

CL0SE(UNIT=3) 

100  F0RMAT(3X,I5,3(3X,F10. 1) ,2(3X,I2) ,3X,F8. 2) 


RETURN 

END 


SUBROUTINE  PAIR2(CRT,KBD,TEMP3,TEMP4) 

**V(r')V***')V**'>Sr***5V****5V**5V*>V**5VA***************yr******i«r*****T>r*Vri>r* 

***  PROGRAM  TO  PAIR  POINT  COUNTS  AFTER  THE  DATA  HAS  BEEN  *** 
***  TESTED  BY  "  RANGE  FOURTH  PASS.  *** 

******V?*****Vf******Vf****************************************** 


DIMENSION  X(200),  Y(200),  Z(200),  SRDIST(200) 
INTEGER*4  PC(200),  ARRAY(200),  MODE(200),  CRT,  HOLD 
CHARACTER*9  TEMP3,  TEMP4 


OPEN( 1,F I LE=TEMP3, STATUS* 'OLD' ) 

OPEN( 2 , FI LE=TEMP4 , STATUS* ' NEW ’ ) 

C 

HOLD  *  0 
IREC  =  0 
NREC  =0 
I  =  1 

C  ...  READ  RECORDS  BY  TWO'S  TO  COMPARE  POINT  COUNTS; 

READ(1,*,END=40,ERR=30)  PC(I),  X(I),  Y(I),  Z(I), 

+  ARRAY(I),  MODE(I),  SRDIST(I) 
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READ(1,*,END=40,ERR=30)  PC(I+1),  X(I+1),  Y(I+1),  Z(I+1), 

+  ARRAY(I+1),  M0DE(I+1),  SRDIST(I+1) 

NREC  =  NREC  +  1 
C 

IF(PC(I+1)  .EQ.  HOLD)  THEN 

WRITE(2,100)  PC(I+1),  X(I+1),  Y(I+1),  Z(I+1), 

+  ARRAY(I+1),  M0DE(I+1),  SRDIST(I+1) 

HOLD  =  PC(I+1) 

GO  TO  20 
END  IF 
C 

IF(PC(I)  .EQ.  PC(I+1))  THEN 

WRITE(2,100)  PC(I),  X(I),  Y(I),  Z(I),  ARRAY(I),  MODE(I), 
+  SRDIST(I) 

WRITE(2,100)  PC(I+1),  X(I+1),  Y(I+1),  Z(I+1), 

+  ARRAY(I+1),  M0DE(I+1),  SRDIST(I+1) 

HOLD  =  PC(I+1) 

IREC  =  IREC  +  1 
GO  TO  20 
END  IF 
C 

IF(PC(I)  .NE.  PC(I+1))  THEN 
PC(1)  =  PC(I+1) 

X(l)  =  X(I+1) 

Y(l)  =  Y(I+1) 

Z(l)  =  Z(I+1) 

ARRAY(l)  =  ARRAY(I+1) 

MODE(l)  =  M0DE(I+1) 

SRDIST(l)  =  SRDIST(I+1) 

I  =  1 
END  IF 
C 

GOTO  10 

20  1  =  1  +  1 
GO  TO  10 


30  WRITECCRT,*)'  THERE  IS  AN  ERROR  IN  THE  DATA  FILE  IN  REC. ' ,NREC 
WRITE(CRT,*)'  ...  OPERATION  TERMINATING  DUE  TO  ERROR.  ' 

STOP 

40  CONTINUE 


CL0SE(UNIT=1) 

CLOSE(UNIT=2) 

100  F0RMAT(1X,I5,2X,F7. 1,2X,2F9. 1,2X,I2,2X,I2,2X,F8. 2) 

RETURN 

END 
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SUBROUTINE  REC( CRT,KBD,TEMP4) 

**A***'>V**ilc********Ve*******V»*******Vr****************-***** 

***  Program  to  produce  the  final  gating  of  Keyport  *** 
***  hydrophone  array  test  data.  *** 

*****'iV********Vr***********'>V******'************'****i^*'»V*»V'»(r* 


INTEGER  PC,  ARRAY,  PCI,  PC2,  Al,  A2,  ARRAYl,  ARRAY2,  CRT 

INTEGER  PCA(IO),  ARRAYA(IO),  HOLD 

DIMENSION  XA(IO),  YA(IO),  ZA(IO),  SRDISTA(IO) 

CHARACTER  0UTFIL*13,  ANS*1,  TEMP4*9,  TEMP1*9, RANGE*? 


RANGE='NANOOSE' 


99  WRITE(CRT,*)’  What  is  the  name  you  wish  to  give  to  ', 

+  '  your  output  file?  ' 

READ(CRT,'(A)' )  OUTFIL 

WRITE(CRT,*) '  Enter  the  numbers  of  the  arrays  to  be  paired:  ' 
READCKBD,*)  Al,  A2 


OPEN( 1 ,FILE=TEMP4,STATUS='0LD' ) 

OPEN( 2 , F I LE= ' TEMP 1 .  DAT ' , STATUS= ' OLD ’ ) 


5  READ(1,*,END=8,ERR=8)  PC,  X,  Y,  Z,  ARRAY,  MODE,  SRDIST 

IF(( ARRAY  .EQ.  Al).  OR.  (ARRAY  .  EQ.  A2))  THEN 
WRITE(2,100)  PC,  X,  y,  Z,  ARRAY, SRDIST 
END  IF 
GOTO  5 

8  CONTINUE 
CLOSE  (UNIT=1) 

. . .  Routine  to  pair  data  again: 

REWIND  2 

OPEN ( 4 , F I LE= ' TEMP2 . DAT ' , STATUS= ' OLD ' ) 

C 

1  =  1 
IFLAG  =  0 
FIRST  =  1 
HOLD  =  0 
M  =  0 
C 

9  READ(2,100,END=40)  PCA(I),  XA(I),  YA(I),  ZA(I), 

+  ARRAYAC I ) ,  SRDISTAC I ) 

IF( FIRST  .EQ.  1)  THEN 
FIRST  =  0 
HOLD  =  PCA(I) 

END  IF 

IF(PCA(I)  .EQ.  HOLD)  THEN 
1  =  1  +  1 
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...  In  cases  where  there  are  three  or  more  reports  for 
...  a  given  point  count,  segregate  by  comparing  SRDIST 
...  if  there  are  more  than  3  reports,  discard  all. 

IF  (M  .EQ.  3)  THEN 

IFCARRAYACM)  .  EQ.  ARRAyA(M-l))  THEN 
IF  (ARRAYA(M)  .  EQ.  ARRAYA(M-2))  THEN 
GOTO  50 
ELSE 

IF  (ABS(SRDISTA(M-2)-SRDISTA(M))  .  LT. 
ABS(SRDISTA(M-2)-SRDISTA(M-l)))  THEN 
WRITE(4,100)  PCA(M),  XA(M),  YA(M),  2A(M), 
ARRAYA(M),  SRDISTA(M) 

WRITE(4,100)  PCA(M-2),  XA(M-2),  YA(M-2),  ZA(M-2), 
ARRAYA(M-2),  SRDISTA(M>2) 

I FLAG  =  1 
ELSE 

WRITE(4,100)  PCA(M-l),  XA(M-l),  YA(M-l),  ZA(M-l), 
ARRAYA(M-l),  SRDISTA(M-l) 
WRITE(4,100)  PCA(M-2),  XA(M-2),  YA(M-2),  ZA(M-2), 
ARRAYA(M-2),  SRDISTA(M-2) 

I FLAG  =  1 
END  IF 
END  IF 
ELSE 

IF  (ARRAYA(M)  .  EQ.  ARRAYA(M-2))  THEN 
IF  (ABS(SRDISTA(M-l)-SRDISTA(M))  .  LT. 

ABS( SRDISTA(M-l) -SRDISTA(M-2) ) )  THEN 
WRITE(4,100)  PCA(M),  XA(M),  YA(M),  ZA(M), 
ARRAYA(M),  SRDISTA(M) 

WRITE(4,100)  PCA(M-l),  XA(M-l),  YA(M-l),  ZA(M-l), 
ARRAYA(M-l),  SRDISTA(M-l) 

IFLAG  =  1 
ELSE 

WRITE(4,100)  PCA(M-l),  XA(M-l),  YA(M-l),  ZA(M-l), 
ARRAYA(M-l),  SRDISTA(M-l) 
WRITE(4,100)  PCA(M-2),  XA(M-2),  YA(M-2),  ZA(M-2), 
ARRAYA(M-2),  SRDISTA(M-2) 

IFLAG  =  1 
END  IF 
ELSE 

IF  (ABS(SRDISTA(M)-SRDISTA(M-l))  .  LT. 

ABS(SRDISTA(M)-SRDISTA(M-2)))  THEN 
WRITE(4,100)  PCA(M),  XA(M),  YA(M),  ZA(M), 
ARRAYA(M),  SRDISTA(M) 

WRITE(4,100)  PCA(M-l),  XA(M-l),  YA(M-l),  ZA(M-l), 
ARRAYA(M-l),  SRDISTA(M-l) 

IFLAG  =  1 
ELSE 

WRITE(4,100)  PCA(M),  XA(M),  YA(M),  ZA(M), 
ARRAYA(M),  SRDISTA(M) 

WRITE(4,100)  PCA(M-2),  XA(M-2),  YA(M-2),  ZA(M-2), 


on  n  n  o  o  on  on 


+  ARRAYA(M-2),  SRDISTA(M-2) 

IFLAG  =  1 
END  IF 
END  IF 

END  IF 


ELSE 

IF  ((M  .EQ.  2) 
VmiTE(4,100) 

WRITE(4,100) 


END  IF 

END  IF 

50  CONTINUE 


25  CONTINUE 

IFLAG  =  0 
PCA(l)  =  PCA(I) 

XA(1)  =  XA(I) 

YA(1)  =  YA(I) 

ZA(1)  =  ZA(I) 

ARRAYA(l)  =  ARRAYA(I) 
SRDISTA(l)  =  SRDISTA(I) 
I  =  2 
M  =  1 

HOLD  =  PCA(l) 

DO  45  L  =  2,10 
PCA(L)  =  0 
XA(L)  =  0 
YA(L)  =  0 
ZA(L)  =  0 
ARRAYA(L)  =  0 
SRDISTA(L)  =  0 
45  CONTINUE 

END  IF 

GOTO  9 

40  CONTINUE 


.AND.  (ARRAYA(M)  .  NE.  ARRAYA(M-l)))  THEN 
PCA(M),  XA(M),  YA(M),  ZA(M) , ARRAYA(M) , 
SRDISTA(M) 

PCA(M-l),  XA(M-l),  YA(M-l),  ZA(M-l), 
ARRAYA(M-l),  SRDISTA(M-l) 


NREC  =  0 
C 

CLOSE  (UNIT=4) 

C 

OPEN(7,FILE='TEMP2.  DAT' ,STATUS='OLD' ) 

OPEN( 9 ,FILE=OUTFIL, STATUS=’ NEW' ) 

C 

WRITE(9,300)  RANGE,  Al,  A2 
C  ...  Read  in  array  data  in  two  record  pairs: 
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10  READ(7,100,END=60,ERR=70)  PCI,  XI,  Yl,  Zl,  ARRAYl,  SRDISTl 

READ(7,100,END=60,ERR=70)  PC2,  X2,  Y2,  Z2,  ARRAY2,  SRDIST2 


IF( ARRAYl  .EQ.  A1  .AND.  ARRAY2  .  EQ.  A2)  THEN 

...  If  arrays  are  in  specified  order  (e. g.  4,5): 
WRITE(9,200)  PCI,  XI,  Yl,  Zl,  X2,  Y2,  Z2 

END  IF 

IF( ARRAYl  .EQ.  A2  .AND.  ARRAY2  .  EQ.  Al)  THEN 

...  If  arrays  are  in  reverse  order  (eg.  5,4): 
WRITE(9,200)  PCI,  X2,  Y2.  Z2,  XI,  Yl,  Zl 

END  IF 


. . .  Increment  record  counter: 
NREC  =  NREC  +  1 


GOTO  10 


70  WRITE(CRT,*) '  There  is  a  bad  record  in  the  file. ' 
60  CONTINUE 


CL0SE(UNIT=2) 
CL0SE(UNIT=7) 
CLOSE ( UN IT=9) 


100  F0RMAT(I5,3(2X,F8.  1) ,2X,I2,2X,F8.  2) 

200  F0RMAT(2X,I5,1X,6(2X,F11.  D) 

300  F0RMAT(16X,A10,2X,I2,3X,I2) 

WRITE(CRT,*)  '  Do  you  want  to  try  another  array  pair?  (Y/N)' 
READ(KBD, '(A)' )  ANS 

IF(ANS  .EQ.  'Y'  .OR.  ANS  .  EQ.  ’y’)  GO  TO  99 


RETURN 

END 
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APPENDIX  B.  PROGRAM  TOE 


PROGRAM  TOE 


C* 

C*  Progranuner 
C* 

C*  Purpose 
C* 

C* 

C* 

C* 

C  *****************1! 


Sukru  KORLU 

To  calculate  and  eliminate  timing  offset  error  by 
least  squares  in  the  keyport  hydrophone  data. 
Output  of  this  program  can  be  read  by  the  program 
keymainS. 


* 

* 

* 

* 

* 

* 

* 


C 

C 

c 

c 


c 


c 

c 

c 

c 

c 


DIMENSION  DEPTH(53),  VEL(53),  L(53),  G(53) 

INTEGER  ARRl,  ARR2,  NL,  PC,  KBD,  CRT 

REAL*8  X3,  Y3,  VO,  VI,  DEPTH,  VEL,  DL,  L,  G,  LXl,  LX2,  LX3,  LYl, 
+  LY2,  LY3,  PI,  P2.  Al,  A2,  GG,  TTHETA,  XTHETA,  YTHETA, 

+  TVEL,  XVEL,  YVEL 

CHARACTER* 12  FNAMEl,  FNAME2,  FNAME3,  FNAME4,  FNAME5,  FNAME6, 

+  FNAME7,  FNAME8 

CHARACTER*!  ANS 

LOGICAL  EXT,  EXT3,  EXT4,  EXTS,  EXT7 
PARAMETER(KBD=5  ,  CRT=6  ,  NL=53) 


C  Define  input/output  files.  .  . 

C 

WRITE(CRT,*)  '  Please  enter  the  name  of  the  range  configuration' , 
+'  file...:  ' 

READ(KBD,'(A)')  FNAMEl 
INQUIREC  FILE=FNAME 1 , EXIST=EXT) 

IFC.NOT.  EXT)  THEN 

WRITE(CRT,*)  '  The  range  config.  file  does  not  exist.  .  .  ' 

GOTO  99 
END  IF 


WRITE(CRT,*)  '  Please  enter  the  name  of  the  DV  file...: 
READ(KBD,'(A)' )  FNAME6 
INQ*JIRE(FILE=FNAME6  ,EXIST=EXT) 

IF(.NOT.  EXT)  THEN 

WRITE(CRT,*)  '  The  DV  file  does  not  exist...  ' 

GOTO  99 
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END  IF 
C 

WRITE(CRT,*)  '  Please  enter  the  name  you  wish  to  give  to  your' , 
+'  output  file. . . :  ' 

READ(KBD,'(A)')  FNAME8 
I NQU I RE ( F I LE=FN AME  8 , EX I ST=EXT ) 

IF(EXT)  THEN 

WRITE(CRT,*)  '  The  output  file  already  exists. . .  ' 

GOTO  99 
END  IF 
C 

10  WRITE(CRT,*)  '  Please  enter  the  name  of  the  input  file...:  ' 
READ(KBD,'(A)')  FNAME2 
INQUIRE(FILE=FNAME2 ,EXIST=EXT) 

IFC.NOT.  EXT)  THEN 

WRITE(CRT,*)  '  The  input  file  does  not  exist, . .  ' 

GOTO  99 
END  IF 
C 

FNAME3  =  'LOCAL. REG' 

FNAHE4  =  ' TRANS 1. REG' 

FNAME5  =  'TRANS 2. REG' 

FNAME7  =  'THETA. REG' 

I NQU I RE ( F I LE=FNAME  3 , E X I ST=EXT3 ) 

I NQU I RE ( F I LE=FNAME4 , EX I ST=EXT4) 

I NQU I RE ( F I LE=FNAME  5 , EX I ST=EXT5 ) 

INQUIRE(FILE=FNAME7,EXIST=EXT7) 

IF(EXT3  .OR.  F.XT4  .OR.  EXT5  .OR.  EXT7)  THEN 

WRITE(CRT,*)  '  Some  of  the  intermediate  files  alredy  exsist,  .  .  ' 
GOTO  99 
END  IF 


Check  input  file  get  ok.  from  the  user.  .  . 

0PEN( 1,FILE=FNAME2,STATUS='0LD' ) 

READ( 1,100)  ARRl,  ARR2 

WRITE(CRT,*)  '  The  arrays  in  the  file  are:  ',ARR1,ARR2 
WRITE(CRT,*)  '  Is  this  the  file  you  want  ?  (Y/N)  ' 
READ(KBD,'(A)')  ANS 

IF(ANS  .NE.  'Y'  .OR.  ANS  .NE.  'y')  THEN 
TOITE(CRT,*)  '  Redefine  input  data  file.  .  .  ’ 
CLOSE(UNIT=l) 

GOTO  10 
END  IF 

CL0SE(UNIT=1) 


Begin  to  call  subroutines.  . . 


CALL  L0CAL( FNAME 1 , FNAME2 , FNAME3 , ARR 1 , ARR2 , X3 , Y3 , CRT) 
C 
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noon  o  o  non  nn  nn  nnnnn 


C 


CALL  TRANSFRM(FNAME3 ,FNAME4 ,FNAME5 .CRT) 
CALL  LSLINE(FNAME6,V0,V1,NL,CRT) 


OPEN ( 2 , F I LE=FNAME  6 , STATUS= ‘ OLD ' ) 


Read  depth  velocity  profile. . . 

WRITE(CRT,*)  '  Setting  layer  end  points  and  gradient  values. . .  ' 
READ(2,*) 

READ(2,*) 

DL  =  25. ODO 


DO  20  M  =  1  ,  52 
READ(2,200)  DEPTH(M) , VEL(M) 

20  CONTINUE 

CLOSE ( UN IT=2) 

DEPTH(53)  =  DEPTH(52)+DL 
VEL(53)  =  2.  0D0*VEL(52)-VEL(51) 


Setting  layer  end  points  and  gradient  values.  . . 

DO  30  M  =  1  ,  52 
L(M)  =  DEPTH(M)-(DL/2. ODO) 

G(N)  =  VEL(M+1)-VEL(M) 

30  CONTINUE 

GG  =  0. ODO 

DO  40  M  =  51  ,  44  ,  -1 
GG  =  GG+(G(M+1)-G(M)) 

40  CONTINUE 

GG  =  GG/8. ODO 
L(53)  =  L(52)+DL 
G(53)  =  G(52)+GG 

DO  50  M  =  1  ,  53 
G(M)  =  G(M)/DL 
50  CONTINUE 

Read  position  of  the  torpedo  according  to  the  left  and  the  rigth  array 
then  call  subroutine  RTRACE  for  raytracing... 

0PEN( 3 ,FILE=FNAME3 , STATUS=’ OLD' ) 

OPEN(4,FILE=FNAME7 ,STATUS='NEW' ) 

C 

WRITE(CRT,*)  '  Reading  the  position  of  the  torpedo. . .  ’ 

C 

60  READ(3,300,END=80,ERR=70)  PC,LX1,LX2,LX3,LY1,LY2,LY3 
C 

WRITE(CRT,*)  '  Multilayer  ray  tracing...  ' 

C 

PI  =  DSQRT((LX1'''^2)+(LX2**2)) 
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P2  =  DABS(X3)-LX3 
A1  =  0.  ODO 
A2  =  DABS(X3) 

C 

CALL  RTRACE(NL,L,VEL,DEPTH,A1,A2,P1,P2,DL,G,V0,V1,TTHETA,TVEL) 
XTHETA  =  TTHETA 
XVEL  =  TVEL 
C 

PI  =  DSQRT((LY1**2)+(LY2**2)) 

P2  =  DABS(Y3)-LY3 
A1  =  0.  ODO 
A2  =  DABS(Y3) 

C 

CALL  RTRACECNL, L,VEL, DEPTH, A1,A2, PI, P2,DL,G, VO, VI, TTHETA, TVEL) 
YTHETA  =  TTHETA 
YVEL  =  TVEL 
C 

WRITE(4,400)  PC, XTHETA, XVEL, YTHETA, YVEL 
GOTO  60 

70  WRITE(CRT,*)  'There  is  an  error  in  the  local  track  file.  . .  ' 
CL0SE(UNIT=3) 

CL0SE(UNIT=4) 

GOTO  99 

80  CL0SE(UNIT=3) 

CL0SE(UNIT=4) 


CALL  CDELTAC  FNAME2 , FNAME4 , FNAME5 , FNAME  7 , FNAME8 , CRT) 


WRITE(CRT,*)  '  Operation  complete.  Program  terminating.  .  .  ' 


STOP 

99  WRITE(CRT,*)  '  OPERATION  TERMINATING  DUE  TO  THE  ERROR  M!  ' 
STOP 


100  FORMAT(28X,I2,3X,I2) 

200  F0RMAT(8X,F4. 0,1X,F9.4) 
300  FORMAT(I5,1X,6(2X,F10. 2)) 
400  FORMAT(I5,3X,4(F15.  9,2X)) 


END 
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SUBROUTINE  LOCAL( NAME 1 , NAME2 .NAMES , ARRl , ARR2 , X3 , Y3 , CRT) 


C* 

* 

C*  Programmer 

Sukru  KORLU 

* 

C* 

* 

C*  Purpose 

To  calculate  local  track  by  using  global  track  and 

* 

c* 

array  locations. 

* 

C* 

* 

C*  Key 

variables 

* 

C* 

NAMEl 

Name  of  the  input  file  which  carries  the  pair  of 

* 

C* 

global  track. 

* 

C* 

NAME  2 

Name  of  the  range  configuration  file. 

* 

c* 

NAMES 

Name  of  the  output  file  which  carries  the  pair  of 

* 

C* 

local  track. 

* 

C* 

ARRl 

Left  array  number. 

* 

C* 

ARR2 

right  array  number. 

* 

C* 

PC 

Point  count 

* 

c* 

GX1-GX3 

Global  track  from  left  array. 

* 

c* 

GY1-GY3 

Global  track  from  rigth  array. 

* 

c* 

X-Y-Z 

Array  location. 

* 

c* 

ARRNUM 

Array  number. 

* 

c* 

X1-X3 

Left  array  location. 

* 

c* 

Y1-Y3 

Rigth  array  location. 

* 

c* 

LX1-LX3 

Local  track  from  the  left  array. 

* 

c* 

LY1-LY3 

Local  track  from  the  rigth  array. 

* 

c* 

Q1,Q2 

Logical  constants. 

* 

c* 

* 

C 

INTEGER  ARRl,  ARR2,  ARRNUM,  Q,  PC,  CRT 
C 

REAL*8  GXl,  GX2,  GX3,  GYl,  GY2,  GY3,  X,  Y,  Z,  XI,  X2,  X3,  Yl, 
+  Y2,  Y3,  LXl,  LX2,  LX3,  LYl,  LY2,  LY3 

C 

CHARACTER* 12  NAMEl,  NAME2,  NAMES 
C 

C  Read  array  locations. . . . 

C 

0PEN(1,FILE=NAME1,STATUS='0LD' ) 

C 

Q=0 

WRITE(CRT,*)  '  Reading  array  locations. . . 

10  READ(1,100,END=30,ERR=20)  X,Y,Z, ARRNUM 

IF( ARRNUM  .EQ.  ARRl)  THEN 
XI  =  X 
X2  =  Y 
X3  =  Z 
Q  =  Q+1 
END  IF 

IF( ARRNUM  .EQ.  ARR2)  THEN 
Yl  =  X 
Y2  =  Y 
Y3  =  Z 


C 

C 
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Q  =  Q+1 
END  IF 
C 

GOTO  10 

20  WRITE(CRT,*)  '  There  is  an  error  in  the  range  configuration' , 

+  '  file...  ' 

CL0SE(UNIT=1) 

GOTO  99 

30  CL0SE(UNIT=1) 

C 

C  Make  sure  that  both  array  locations  were  found.  . . 

C 

IF(Q  .NE.  2)  THEN 

WRITE(CRT,*)  '  At  least  one  array  location  can  not  be  found' , 
+  '  in  the  range  configuration  file.  . .  ' 

CL0SE(UNIT=1) 

GOTO  99 
END  IF 
C 

C  Calulate  local  track. .  , 

C 

OPENC  2 ,FILE=NAME2 , STATUS=' OLD ' ) 

0PEN( 3 ,FILE=NAME3 , STATUS=' NEW* ) 

C 

WRITE(CRT,*)  '  Reading  global  track  writing  local  track... ' 
READ(2,*) 

C 

40  READ(2,200,END=60,ERR=50)  PC ,GX1 ,GX2 ,GX3 .GYl ,GY2 ,GY3 
LXl  =  GXl-Xl 
LX2  =  GX2-X2 
LX3  =  GX3-X3 
LYl  =  GYl-Yl 
LY2  =  GY2-Y2 
LY3  =  GY3-Y3 

WRITE(3,300)  PC,LX1,LX2,LX3,LY1,LY2,LY3 
GOTO  40 

50  WRITE(CRT,*)  '  There  is  an  error  in  the  input  file.  .  .  ' 
CL0SE(UNIT=2) 

CL0SE(UNIT=3) 

GOTO  99 

60  WRITE(CRT,*)  '  Transformation  done...  ' 

C 

CL0SE(UNIT=2) 

CLOSE ( UN IT=3) 

C 

RETURN 

C 

C 

99  WRITE(CRT,*)  'OPERATION  TERMINATING  DUE  TO  THE  ERROR  !!!  ’ 

STOP 

C 

100  F0RMAT(F10.  4,3X,F11.  3,6X,F11.  3,6X,I2) 

200  F0RMAT(2X,I5,1X,6(2X,F11.  D) 

300  F0RMAT(I5,1X,6(2X,F10. 2)) 

C 
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END 

C 

C 

c 

SUBROUTINE  TRANSFRM( NAME 1 ,NAME2 , NAMES , CRT) 

C'>V*iV******'>V*Vf)V***'>VA’**iV*>V*********Vfilr>V****ilr***-**iV****i>r**********ifr***'****** 


c* 

* 

C*  Programmer 

Sukru  KORLU 

★ 

c* 

* 

C*  Purpose 

To  make  transform  from  cartesian  coordinates  to 

the  * 

c* 

azimuth, horizontal  and  vertical  components. 

* 

c* 

* 

C* 

Key  Variables 

* 

C* 

NAMEl 

Name  of  the  file  that  carries  the  pair  of  local 

track. * 

c* 

NAME2-3 

Names  of  the  output  files. 

it 

c* 

PC 

Point  count. 

* 

c* 

LX1-LX3 

Local  track  from  the  left  array. 

* 

C* 

LY1-LY3 

Local  track  from  the  rigth  array. 

* 

c* 

FIX 

Azimuth  component  of  the  left  array  track. 

* 

c* 

FIY 

Azimuth  component  of  the  rigth  array  track. 

it 

c* 

HORX 

Horizontal  component  of  the  left  array  track. 

* 

c* 

HORY 

Horizontal  component  of  the  rigth  array  track. 

* 

c* 

VERX 

Vertical  component  of  the  left  array  track. 

it 

c* 

VERY 

Vertical  component  of  the  rigth  array  track. 

it 

c* 

* 

C*A*****,V**************************************************************Vc 

c 

INTEGER  PC,  CRT 
C 

REAL*8  LXl,  LX2,  LX3,  LYl,  Ly2,  LY3,  FIX,  HORX,  VERX,  FIY, 

+  HORY,  VERY,  PI 

C 

CHARACTER*12  NAMEl,  NAME2,  NAMES 
C 

PARAMETER(PI=3.  141592654D0) 

C 

OPEN( 1 , FILE=NAME 1 , STATUS= ' OLD ' ) 

OPEN( 2 , FILE=NAME2 , STATUS=' NEW' ) 

OPEN( 3 ,FILE=NAME3 , STATUS=’ NEW’ ) 

C 

WRITE(CRT,*)  '  Making  tranform.  .  .  ’ 

C 

10  READ(1,100,END=30,ERR=20)  PC,LX1,LX2,LX3,LY1,LY2,LY3 
FIX  =  DATAN(LX2/LX1) 

IF(LX1  .LT.  0)  FIX=FIX+PI 
HORX  =  LX1/DC0S(FIX) 

VERX  =  LX3 

FIY  =  DATAN(LY2/LY1) 

IF(LY1  .LT.  0)  FIY=FIY+PI 
HORY  =  LYl/DCOSCFIY) 

VERY  =  LY3 

WRITE(2,200)  PC, FIX, HORX, VERX 
WRITE(3,200)  PC, FIY, HORY, VERY 
GOTO  10 

20  WRITECCRT,*)  '  There  is  an  error  in  the  local  track  file...' 
CLOSE(UNIT=l) 
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CL0SE(UNn=2) 

CL0SE(UNIT=3) 

GOTO  99 

30  CL0SE(UNIT=1) 

CL0SE(UNIT=2) 

CL0SE(UNIT=3) 

C 

RETURN 

C 

C 

99  WRITE(CRT,*)  ’OPERATION  TERMINATING  DUE  TO  THE  ERROR  !!!  ' 
STOP 
C 

100  F0RMAT(I5,1X,6(2X,F10. 2)) 

200  FORMAT(I5,3X,3(F15.  8,2X)) 

C 

C 

END 


C 

C 

C 


SUBROUTINE  LSLINE( FNAME , VO , VI , NL , CRT) 


Programmer 


C* 

C* 

C* 

C*  Purpose 
C* 

C* 

C*  Key  variables 


Sukru  KORLU 

To  calculate  VO  (speed  of  sound  at  the  surface)  and 
VI  (gradient)  by  least  squares. 


* 

C* 

FNAME 

C* 

c* 

X 

c* 

SX 

c* 

NL 

c* 

XB 

c* 

SSX 

c* 

Y 

c* 

SY 

c* 

YB 

c* 

SXY 

c* 

VO 

c* 

VI 

c* 

r***y 

C 

C 

C 


Name  of  the  file  which  carries  depth  velocity 
profile. 

Depth. 

Sum  of  X. 

#  of  layers. 

Mean  of  x. 

Sum  of  square  of  x. 

Speed  of  sound  in  the  water. 

Sum  of  y. 

Mean  of  y. 

Sum  of  product  of  x  and  y. 

Sound  velocity  at  the  surface. 

Gradient. 


INTEGER  NL,  CRT 

REAL*8  VO,  VI,  X,  SX,  SSX,  XB,  Y,  SY,  YB,  XY,  SXY 
CHARACTER  FNAME* 12 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


C 

0PEN( 1 , FILE=FNAME , STATUS=' OLD’ ) 

C 

WRITE(CRT,*)  ’  Calculating  VO, VI  ...  ’ 
READd,*) 
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READd,*) 

C 

SX  =  0.  ODO 
SSX  =  0.  ODO 
SY  =  0.  ODO 
SXY  =  0.  ODO 
C 

DO  10  I  =  1  ,  NL-1 
READd,  100)  X.Y 
SX  =  SX+X 
SSX  =  SSX+(X**2) 

SY  =  SY+Y 
SXY  =  SXY+X*Y 
10  CONTINUE 
C 

XB  =  SX/(NL-1) 

YB  =  Sy/(NL-1) 

VI  =  (SXY-((NL-1)*XB*YB))/(SSX-((NL-1)*(XB**2))) 
VO  =  YB-(V1*XB) 

C 

CL0SE(UNIT=1) 

C 

RETURN 

C 

100  F0RMAT(8X,F4.  0,1X.F9,4) 

C 

END 


C 

C 

c 

c 


SUBROUTINE  RTRACE(NL,L,V,DEPTH,A1,A2,P1,P2,DL,G,V0,V1,TTHETA,TVEL) 

C***************Vr**Vr'jWnV**Vf******Vf******yr******iV****ilf*****ilr************** 

C*  * 

C*  Programmer  :  Sukru  KORLU  * 


C* 


* 


c* 

Purpose  : 

To  calculate  elevation  angle  at  the  torpedo 

by 

* 

c* 

multilater  raytracing  for  one  point  count. 

* 

c* 

Also  to  find  speed  of  sound  at  the  depth  of 

* 

c* 

torpedo. 

* 

c* 

* 

c* 

Key  variables 

* 

c* 

L  : 

An  array  which  indicates  layer  end  points. 

* 

c* 

V  : 

An  array  which  indicates  speed  of  so\ind. 

* 

c* 

G  : 

An  array  which  indicates  gradient  values. 

* 

c* 

THETA  : 

An  array  which  indicates  elevation  angle  in 

each 

* 

c* 

layer. 

* 

c* 

RV  : 

Ray  invariant. 

* 

c* 

R  : 

Range  of  torpedo  from  the  array. 

* 

c* 

A1,A2 

Array  position. 

* 

c* 

P1,P2  : 

Torpedo  position. 

★ 

c* 

VO, VI  : 

Constants. 

* 

c* 

C1,C2  : 

Center  of  the  circle  that  we  assumed  sound  travels 

* 

c* 

on  it's  arc. 

* 

c* 

* 
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DIMENSION  LA(53),  L(53),  VE(53),  V(53),  DEPTH(53),  G(53),  GR(53), 
+  THETA(53),  VVO(53),  C2(53),  T(53) 

C 

INTEGER  NL 
C 

REAL*8  EP,  LA,  L,  VE,  V,  Al,  A2,  PI,  P2,  G,  VO,  VI,  GR,  TVEL,  Cl, 
+  C2,  CC2,  THETA,  ATHETA,  TTHETA,  OTHETA,  WO,  T,  RV,  R,  DL, 

+  DEPTH 

C 

C  Find  the  layers  that  torpedo  or  array  in  it,  redefine  the  end  points 
C  of  those  layers,  define  tolerance... 

EP  =  0. lD-6 
1  =  0 
J  =  0 
C 

DO  10  K  =  1  ,  NL 
LACK)  =  L(K) 

VE(K)  =  V(K) 

GR(K)  =  G(K) 

IFCLACK)  .LE.  A2)  J=J+1 
IF(LA(K)  .LE.  P2)  1=1+1 
10  CONTINUE 
C 

IFCDEPTHCI)  .LE.  P2)  THEN 

VE(I)=VE(I)+GR(I)*(P2-DEPTH(I)) 

ELSE 

VE( I)=VE( I ) -GR( I - 1 )*( DEPTHC I ) -P2 ) 

END  IF 

j  C 

IFCDEPTHCJ)  .JE.  P2)  THEN 

VE(J)=VE(J)+GR(J)*(A2‘DEPTH(J)) 

ELSE 

VE(J)=VE(J)-GR(J-1)*(DEPTH(J)“A2) 

END  IF 
C 

LAC  I)  =  P2 
LACJ)  =  A2 
TVEL  =  VECI) 

C 

C  Calculate  an  initial  estimate  for  the  elevation  angle  using  single 
C  layer  approximation. . . 

CC2  =  -VO/Vl 

Cl  =  CO.  5D0*CP1+A1))+C0.  5D0*C(LACI)-LACJ))*C(CLACI)+LACJ)-2.  ODO* 
+CC2))/CP1-A1)))) 

ATHETA  =  DATANCCA1"C1)/CA2-CC2)) 

C 

C  Using  this  initial  estimate  of  the  elevation  angle  with  ray  invariant 
C  to  raytrace  back  through  all  the  layers... 

20  RV  =  DCOSC ATHETA) /VECJ) 

DO  30  K  =  I  ,  J 
THETACK)  =  DAC0SCRV*VECK)) 

TCK)  =  DTANCTHETACK)) 

VVOCK)  =  VECK)-LACK)*GRCK) 

C2CK)  =  -W0CK)/GRCK) 
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30  CONTINUE 


C 

C 

C 


Using  the  angle  just  calculated  iterate  backwards  through  the  layers 
from  array  to  torpedo  to  get  the  horizontal  range. .  . 

R  =  0.  ODO 

DO  40  K  =  J  ,(I+1),  -1 
Cl  =  R-T(K)*(LA(K)-C2(K-1)) 

R  =  C1+T(K-1)*(LA(K-1)-C2(K-1)) 

40  CONTINUE 

Test  if  the  value  for  the  range  within  tolerance  if  not, redefine  ATHETA 
the  initial  angle  and  raytrace  again. . . 

IF((DABS(R-P1))  .LE.  EP)  GOTO  50 
OTHETA  =  ATHETA 

ATHETA  =  DATAN(DTAN(0THETA)*(R/P1)) 

GOTO  20 


50  TTHETA  =  THETACl) 
RETURN 


C 

C 

c 


END 


SUBROUTINE  CDELTA( NAME 1 , NAME2 , NAME3 , NAME4 , NAME5 , CRT) 


CiVyr/rAAiViV*********-!’?****************************************************** 

C* 

* 

C*  Programmer 
r* 

Sukru  KORLU 

* 

* 

C*  Purpose 

To  calculate  timing  offset  error( delta)  by 

least 

* 

c* 

squares. 

* 

c* 

* 

C*  Key 

variables 

* 

C* 

NAMEl 

Name  of  the  file  which  carries  global  track 

from 

* 

c* 

pair  of  arrays. 

* 

c* 

NAME2 

Name  of  the  file  which  carries  transformed 

track 

* 

c* 

from  left  array. 

* 

c* 

NAME  3 

Name  of  the  file  which  carries  transformed 

track 

* 

C* 

from  rigth  array. 

* 

c* 

NAME4 

Name  of  the  file  which  carries  elevation  angle  and 

* 

c* 

sound  velocity  records. 

* 

c* 

NAME5 

Output  file  which  carries  corrected  track. 

* 

c* 

SSQ 

Sum  of  squares. 

* 

c* 

SPR 

Sum  of  products. 

* 

c* 

X 

Track  from  the  left  array. 

* 

c* 

Y 

Track  from  the  rigth  array. 

* 

c* 

FIX-FIY 

Azimuth  angles. 

* 

c* 

XVEL-YVEL 

Speed  of  sound  at  the  depth  of  torpedo. 

* 

c* 

X/YTHETA 

Elevation  angles. 

* 

c* 

DELTA 

Timing  offset  error. 

* 

c* 

* 

C 

C 
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no  non 


DIMENSION  X(3),  Y(3),  XY(3),  A(100,3),  B(100,3),  AB(3) 

C 

INTEGER  PC,  CRT 
C 

REAL*8  X,  Y,  FIX,  FIY,  XTHETA,  YTHETA,  XVEL,  YVEL,  A,  B,  AB, 
+  XY,  SSQ,  SPR,  DELTA 

C 

CHARACTER*12  NAMEl,  NAME2,  NAME3,  NAME4,  NAMES 
C 

CHARACTER*35  AUTL 


WRITECCRT,*)  '  Calculating  delta. . .  ’ 

OPENC 1 ,FILE=NAME1 , STATUS=’ OLD' ) 

OPEN( 2 ,FILE=NAME2 , STATUS=' OLD' ) 

OPENC  3 ,FILE=NAME3 , STATUS=' OLD' ) 
0PEN(4,FILE=NAME4,STATUS='0LD' ) 


I  =  0 

SSQ  =  0.  ODO 
SPR  =  0.  ODO 
C 

READCl,*) 

C 

10  READCl, 100, END=30,ERR=20)  XC 1) ,XC2) ,XC 3) ,YC 1) ,YC 2) ,YC 3) 
I  =  I+l 

READC2,200)  FIX 
READC3,200)  FIY 

READC4,300)  XTHETA , XVEL , YTHETA , YVEL 

XYCl)  =  XCD-YCl) 

XYC2)  =  XC2)-YC2) 

XYC3)  =  XC3)-YC3) 

AC  1,1)  =  XVEL*DCOSCXTHETA)*DCOSCFIX) 

AC  1,2)  =  XVEL*DCOSCXTHETA)*DSINCFIX) 

AC  1, 3)  =  XVEL*DSINC XTHETA) 

BCI,1)  =  YVEL*DCOSCYTHETA)*DCOSCFIY) 

BCI,2)  =  YVEL*DCOSCYTHETA)*DSINCFIY) 

BCI,3)  =  YVEL*DS INC YTHETA) 

ABCl)  =  ACI,1)-BCI,1) 

ABC2)  =  ACI,2)-BCI,2) 

ABC3)  =  ACI,3)-BCI,3) 

SSQ  =  SSQ+CABC1)**2)+CABC2)**2)+CABC3)**2) 

SPR  =  SPR+CXYCl)*ABCl))+CXyC2)*ABC2))+CXYC3)*ABC3)) 

GOTO  10 

20  WRITECCRT,*)  '  There  is  an  error  in  the  input  file. . .  ' 
CLOSE  CUNIT=1) 

CLOSE  CUNIT=2) 

CLOSE  CUNIT=3) 

CLOSE  CUNIT=4) 

GOTO  99 

30  CL0SECUNIT=2) 

CLOSECUNIT=3) 

CL0SECUNIT=4) 

C 
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DELTA  =  SPR/SSQ 

WRITE(CRT,*)  ’  TIMING  ERROR  =  ' .DELTA 
C 

REWIND  1 
C 

OPEN( 7 ,FILE=NAME5 ,STATUS='NEW' ) 

C 

WRITE(CRT,*)  '  Writing  output  file. . .  ' 

C 

READ( 1,400)  AUTL 
WRITE(7,400)  AUTL, DELTA 
C 

J  =  0 

40  READ(1,500,END=60,ERR=50)  PC,X( 1) ,X(2) ,X(3) ,Y( 1) ,Y(2) ,Y(3) 

J  =  J+1 

X(l)  =  X(1)-DELTA*A(J,1) 

X(2)  =  X(2)-DELTA*A(J,2) 

X(3)  =  X(3)-DELTA*A(J,3) 

Y(l)  =  Y(1)-DELTA*B(J,1) 

Y(2)  =  Y(2)-DELTA*B(J,2) 

Y(3)  =  Y(3)-DELTA*B(J,3) 

WRITE( 7,500)  PC,X(1).X(2),X(3),Y(1),Y(2),Y(3) 

GOTO  40 

50  WRITE(CRT,*)  '  There  is  an  error  in  the  input  file. .  .  ' 
CLOSE  (UNIT=1) 

CLOSE  (UNIT=7) 

GOTO  99 

60  CL0SE(UNIT=1) 

CL0SE(UNIT=7) 

C 

RETURN 

C 

99  WRITECCRT,*)  '  OPERATION  TERMINATING  DUE  TO  THE  ERROR  !!!  ' 
STOP 
C 

100  F0RMAT(8X,6(2X,F11.  D) 

200  FORMAT(8X,F15.  8) 

300  F0RMAT(8X,4(F15.  9,2X)) 

400  FORMAT(A35,5X,F8.  6) 

500  F0RMAT(2X,I5,1X,6(2X,F11.  1)) 

C 

END 
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APPENDIX  C.  APL  CODE  TO  CALCULATE  ARRAY  DISPLACEMENT 
AND  ORIENTATION  CORRECTIONS 


-  V-LOC  CORE  SET 
Cl]  EPSA*0. 01 
C2]  X^SETL;  234] 

C3]  Y<-SETLt  5  6  7] 

[4]  XB^  3  1  p((+/X')*N^l^pX) 

[5]  YB^  3  1  p(<+/y)+iV) 

[6]  X-«-X-(pX)pXB 

[7]  y**-y-(py)pyB 

[8]  CXX-f-((W)+.xX)*i7-l 
C9]  cyx<-((isjy)+.xx)tiv-i 
CIO]  cyy-*-((isiy)+.xy)*jv-i 

Cll]  3  3  p  1  0  0  0 

Cl 2]  D^YB-XB 

C13]  BC3;l]-«-0 

C14]  A*-D-iB-I)  +  .xLOC 

C15]  (+/((y-X+.xW)*2))*iy 

C16]  (yB-A+B+.xXB)*2 

C17]  IQ^+/IQ0*1QL 

C18]  L:Zt-CYX+iiYB~A)+.x^XB') 

C19]  BB*B  BMAXl  Z 
C20]  Bt-BB 

C21]  BB-«-yB-(ZiOC+B+.x(XB-LOC)) 

C22]  BBC3;l]-»-0 

C23]  BA+(+/(U-AA)*2))*0.5 

C24]  D*DD 

C25]  *Lii\EPSA<\RA 
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C263  A^D~(.B-I)^.'>^LOC 

[27]  fll>*-3pO 

[28]  LCO-*-,  (+/((Y-X+.xlsjB)*2))*iV 

[29]  LCL-*-,  (yB-A+B+.xXB)*2 

[30]  Bl[2]-«-“loB[3;l] 

[31]  Bl[l]-«-”lo(B[3;2]*(2oBl[2]  )) 

[32]  Bl[3]-^“lo(B[2;l]*<2oBl[2] )) 

[33]  A-*r,A 

[34]  LQ^-^/LQO^LQL 

[35]  7-*-,B,Bl 
V 


-  BB-B  BMAX  D 

[1]  EPSB*1E~B 

[2]  W“*-V0++/  +  /Dx^B 

[3]  LiWW^W 

[4]  C-*-D+.^^B 

[5]  BB^TWODIM  C 

[6]  BB-*-BB+.xB 

[7]  W^+/+/Dx^BB 

[8]  R-*-W-WW 

[9]  ^Lxi£'PSB<|i? 

7 
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V  BB^TWODIM  D;SiCiSSiBltEiB2iD 

111  S**-Z)C3;2]-DC2;3] 

C2]  C>*-£JC2;2]+£)[3;3] 

C3]  5S-«-((5xS)+C5<C)*0.5 

[4]  S^S*SS 
C5]  C<rC*SS 

[6]  Bl-*-  3  3  p  1  0  0  0  ,C,S,0,  <-S),( 

[7]  E-«-D+.xB1 

C8]  S>*-EC3;1]-EC1;3] 

C93  C**-E[3;3]+£’C1;1] 

[10]  5S'«-(  (Sx5)+CxC)*0. 5 

Cll]  S^SiSS 

[12]  C^C*SS 

[13]  B2-«-  3  3  PC,0,5,0,1,0,(-S),0,C 

[14]  f-«-£:+.xB2 

[15]  S-^F[2  5l]-F[l;2] 

[16]  C-«-F[l;l]+F[2;2] 

[17]  SS-f-C  (5xS)+CxC)*0. 5 

[18]  S^S-iSS 

[19]  C'rCiSS 

[20]  B3'e-  3  3  pC,S ,  0  ,  ( -S  )  ,C,  0  0  0  1 

[21]  BB'«-6(B1  +  .xS2  +  .xB3 

V 


F;B3 
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APPENDIX  D.  APL  CODE  TO  CALCULATE  NONNEGATIVE 
COVARIANCE  MATRICES 

-  IND  REGRESS  DATA 
Cl]  INDl^WIND 

[2]  X<ri+/DATA:*(l‘^pDATA) 

[3]  Y**-((plJVi31),6)pO 

[4]  MSE^MSB^t-  6  6  pO 

C5]  D*-0 

C6]  Jfl 

C7]  LliT^D+R-t-INDlLJl 

C8]  YC J;  ]’*-(+/(  (D,0)  +  ((Y. 6  )))♦;? 

C9]  MSB-*-MSB+((YCJ;]-Y)o.x(YCJ;]-Y)) 

CIO] 

Cll]  i2:Yl-«-ByirACI;]-YCJ{] 

C12]  MSE^MSE+aio.xYl) 

C13] 

C14]  *U^T)/L2 

Cl  5]  D^T 

C16]  J^J+1 

C17]  ->-(.J£pINDl)/Ll 

CIS]  MSB^MSB*DFB^(.pJNDl)~l 

C19]  MSE^MSE*DFE*(.+/INDl)-(.{ri^pIND'>-l'>*{il^pIND')-l'>) 

C20]  C^(  (+/IiVBl) -((  +  /(+/ (IiVB*2  )  )  )*i+/INDl  )  )  )♦  (  (pJ^VBl  )-l  ) 

C21]  EIElt-i(.,EIEt-  1  6  ^EEtrEIGENR  MSB)+1B"20  )*"0 . 5 

-22-  DE-  6  6  -BIB1-1-, 0, 0, 0, 0,0,0,EIBl-2-, 0,0, 0,0, 0,0,EJBl-3-, 0,0,0, 
0,0,0,BJBl-4-,0,0,0,0,0,0,BIBl-5-,0,0.0,0,0,0,BJEl-6- 
C23]  QE^  1  0  +BF 
C24]  LE-*-(QE+.xPE)  +  .x(^(JE) 

C25]  L-*-(<5BB)  +  .x(MSB+.xIiB) 

C26]  EIL^Al  6  ^EL^EIGENR  L) 
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[26]  6  ^EL^EIGENR  L) 

[27]  K^+/A^EIL^1 

[28]  QL*-  1  0  *EL 

[29]  P^(^(fSiLE'>)  +  ,xQL 

[30]  EILl<rA/EIL 

[31]  -»-(X>0)/L3 

[32]  SIGMAB^  6  6  pO 

[33]  -»-L4 

[34]  L3:PK^A/P 

[3  5]  LAMDA<-IDENT-^(K,K)pltB-t-KpO 

[36]  X<-1 

[37]  LSlLAMDALXil^IDENTLXilxEILllXl 

[3  8]  X<*-X+l 

[39]  ^(X^K)/L5 

[40]  SIGMABt-  iPK+.x((  LAMDA -I  DENT  )  + .  x  (  ^PK  )))*C 

[41]  L^iSIGMAE^i(DFBx(MSB-(CxSIGMAB)) )+(DFExMSE))*(DFB+DFE) 

V 
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APPENDIX  E.  APL  CODE  TO  ESTIMATE  STABLE  CORRECTIONS 

7  ALPHA^EST 

-1-  R  Y  .SIGHAB  »SIGVIAE  ARE  GLOBAL  VECTORS 

[2]  Xl-*-  7  2  p  1  0 

[3]  X2-«-  2  2  p  0  1 

[4]  X3-e  1  2  p  1  0 

C5]  Jl-*-  7  7  pi 

[6]  J2-*-  2  2  pi 

[7]  J3-*-  1  1  pi 

C8]  A-*rPERROSE  SIGVIAE 

C  9 ]  C1+ (PENROSE (SIGMAE+ ( 7  xSIGMAB ) ) )  + . xSIGMAB* . 

CIO]  Z?2-^  (PFiV^OSE  (SICM>1£+  (  2  XSIGMAB  )  )  )  + .  xSIGMAB* .  xA 
C 1 1  ]  D3>«-  (PENROSE  (SIGMAE*SIGMAB  )  )  + .  xSICMAB+ .  x;i 

C12]  Z-t-  6  6  pO 
C13]  XlA^iSf((>S)A'>,Z'> 

C14]  Xlil-«-XlA,XlA,Xl>l,XlA,XlA,XlA,XlA 
CIS]  X2A*<H(Z,(^A)) 

C16]  X2i5-*-X2>l,X2;i 

C17]  X3A*«-<S}((W),Z) 

C18]  X,71-*-(«5Xl  )  +  .x,71 

C19]  XJ2^<W2)  +  .xJ2 

C20]  XJ3-«-(«5X3)  +  .xJ3 

C21]  XJDl^fsi((inDl),Z) 

C22]  X,7iJl-»-X.7Pl  ,XJD1  ,XJD1  ,XJD1  ,XJD1  ,XJD1  ,XJD1 
C23]  XJD2^fSt(Z,(^D2)) 

C24]  XJD2^XJD2 »XJD2 
C25]  X.7P3«*-^((HP3),Z) 

C26]  T2‘f-(X1A-XJD1)AX2A-XJD2),(X3A-XJD3'> 

C27]  XXl-«-(<SjXl)  +  .xXl 
C28]  XX2*(l5iX2)  +  .xX2 
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C29]  XX3*e(«»X3)  +  .xX3 

C30]  XX1A-^(I?((I5(7xA)),Z)),(«S(Z.Z)) 

[31]  XX2A-*-(lSj(Z,Z)),<6»(Z,(N<2Xil)))) 

[32]  XX3>l*f(l?((W),Z)),  (<9(Z.Z)) 

[33]  XJXl-*-(Wl)  +  .xJl  +  .xXl 

[34]  XJX2**-(W2)  +  .xJ2  +  .xX2 

[35]  XJX3-«-(W3  )  +  .xJ3  +  .xX3 

[36]  XJXZ)l-^(^((^(49xDl)),Z)),(^(Z,Z)) 

[37]  XJXD2^(lSl(Z,Z)),(l5i(Z,(f5i(HxD2)))) 

[38]  XJXZ)3-^(S>((SJD3),Z)),(«?(Z.Z)) 

[39]  Z’l-*-<3x(XXlil+XX2A+XX3A))-(XJXDl+XJXD2  +X JXD 3  ) 

[40]  Z’ll>«-ri[i6;  i6] 

[41]  ri2*«-ri[6  +  i6  j6  +  a6] 

[42]  2’ll[l;3]-*-2’12[l;3]t-0 

[43]  ITll^PENROSE  Til 

[44]  IT12fPENROSE  T12 

[45]  ALPffA-«-(I2’l^(^((l5l2’ll),Z)),(^(Z.(l5J2’12))))  +  .xr2  +  .xY 
7 
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